NASA Contractor Report 198416 


// ^ 


/ / o 


f 


n 

- 


Numerical Prediction of Turbulent Oscillating 
How and Heat Transfer in Pipes With 
Various End Geometries 

(NASA-CR-198416) NUMERICAL N96-18503 

PREDICTION OF TURBULENT OSCILLATING 
FLOW AND HEAT TRANSFER IN PIPES 

WITH VARIOUS END GEOMETRIES Ph.O. Unclas 

Thesis, Final Report (Minnesota 
Uni v. ) 230 p 

G3/34 0098928 


Kirk LeRoi Oseid 
University of Minnesota 
Minneapolis, Minnesota 


December 1995 


Prepared for 

Lewis Research Center 

Under Grant NAG3-1024 



National Aeronautics and 
Space Administration 


NUMERICAL PREDICTION OF TURBULENT OSCILLATING FLOW AND 
HEAT TRANSFER IN PIPES WITH VARIOUS END GEOMETRIES 


A Thesis 

Submitted to the Faculty of the Graduate School 
of the University of Minnesota 

by 


Kirk LeRoi Oseid 


in Partial Fulfillment of the Requirements 
for the Degree of 
Doctor of Philosophy 


October 1993 



ACKNOWLEDGMENTS 


Producing this thesis has given me a great deal of satisfaction as it marks the end 
of a long, frustrating and in the end very satisfying effort. Four years ago I had only 
a basic understanding of numerical methods. While my interest in Computational 
Fluid Dynamics and heat transfer had peaked, I had little opportunity to perform 
this type of work on a productive and professional level. It became obvious that 
formal graduate studies were required to provide me with the credibility needed to 
practice in this field. I would like to thank my advisor, Dr. Suhas V. Patankar, for 
providing the guidance necessary to study such complex topics. Perhaps even greater 
thanks are in order to Dr. Patankar for taking me on as an advisee in the first place. 
I worked for some months before fully appreciating this facet of my good fortune. 

I also thank Dr. Terry Simon, director of the experimental oscillating flow re- 
search at the University of Minnesota, for his encouragement and interest in my 
work, and for providing access to his own oscillating flow and heat transfer library. 

Dr. James Hodges, whose fascination with fluid dynamics and heat transfer was 
infectious, is responsible for my own interest in these topics. 

I am indebted to John Chai for the generous assistance he offered me during the 
course of my degree program at the University of Minnesota. 

The facilities and computing resources of the Minnesota Supercomputer Institute 
are gratefully acknowledged. I hope always to have similarly comfortable surround- 
ings and modem equipment available to me in the future. 

Special thanks are in order to my parents and parents-in-law for their enthusiastic 
support during my studies. Their faith in me has far surpassed my own. 

My greatest thanks are reserved for my wife Mary. She shared my commitment 
to return to school and cheerfully accepted the role of provider. Her love and en- 
couragement are present in this dissertation. 



Contents 


1 Introduction 1 

1.1 Invention of the Stirling Engine 1 

1.2 Renaissance of the Stirling Engine 3 

1.3 The Stirling Engine in Space 5 

1.4 Other Stirling Engine Applications 6 

1.5 Literature Survey of Oscillating Flow and Heat Transfer 8 

1.5.1 Laminar Oscillating Flow 8 

1.5. 1.1 Experiments 8 

1.5. 1.2 Analytical and Numerical Solutions 9 

1.5.2 Transitional and Turbulent Oscillating Flow 9 

1.5.2. 1 Experiments 19 

1. 5.2.2 Numerical Solutions 12 

1.5.2. 3 Turbulence Modeling 13 

1.5.3 Heat Transfer 17 

1.5. 3.1 Experiments 17 

1.5.3. 2 Analytic and Numerical Solutions 19 

2 The Mathematical Representation of the Problem 21 

2.1 Laminar Flow 21 

2.2 Turbulent Flow 22 

i 



2.2.1 Turbulence Modeling 23 

2. 2. 1.1 The Mixing- Length Model 23 

2. 2. 1.2 The One- Equation Models 24 

2.2. 1.3 The Two-Equation Models 25 

2. 2. 1.4 A Multiple Time Scale Model 30 

2. 2. 1.5 A Renormalization Group Model 31 

2.2. 1.6 Other Turbulence Models 33 

2.2.2 The Mean Flow 35 

3 The Laminar Flow Solution 38 

3.1 Fully Developed Pipe Flow 38 

3.2 Oscillating Pipe Flow 40 

3.2.1 Analytic Solution 40 

3.2.2 Numerical Solution 41 

4 The Turbulent Flow Solution 58 

4.1 Fully Developed Pipe Flow 59 

4.2 Oscillating Pipe Flow 61 

4.2.1 Results for the Sudden Expansion/ Contraction Ends .... 62 

4.2.2 Results for Smooth Nozzle Ends 67 

4.2.2. 1 Shortcomings of the Lam-Bremhorst Model .... 69 

4.2.3 Results of the Kim Model 70 

4.2.4 Results of the RNG Model 74 

5 The Laminar Heat Transfer Solution 128 

5.1 Fully Developed Laminar Heat Transfer 128 

5.2 Heat Transfer in Laminar Oscillating Pipe Flow 130 


u 



6 The Turbulent Heat Transfer Solution 163 

6.1 Fully Developed Turbulent Heat Transfer 163 

6.2 Heat Transfer in Turbulent Oscillating Pipe Flow 164 

7 Concluding Remarks 202 

7.1 Summary and Conclusions 202 

7.2 Contributions of This Work 203 

7.3 Suggestions for Further Research 204 


in 



List of Figures 


2. 1 Axisymmetric domain geometry and computational grid for pipe with 
expansion/contraction ends 


37 


3.1 Fully developed velocity profiles given by the numerical solution for 

three grid densities 

3.2 Velocity profiles throughout the oscillating flow cycle for Va = 1 as 

given by Sexl’s analytic solution 

3.3 Velocity profiles throughout the oscillating flow cycle for V a = 30 as 

given by Sexl’s analytic solution 

3.4 Velocity profiles throughout the oscillating flow cycle for Va = 100 

as given by Sexl’s analytic solution 

3.5 Variation of the bulk velocity-pressure gradient phase shift with Va. 

3.6 Axisymmetric domain geometry and computational grid for pipe with 

expansion/contraction ends 

3.7 Comparison of exact and finite volume solution velocity profiles for 

Va = 30 

3.8 Velocity profiles at the pipe midlength at select times during the flow 

cycle for Va = l 

3.9 Velocity profiles at the pipe midlength at select times during the flow 

cycle for Va = 30 


44 

45 

46 

47 

48 

49 

50 

51 

52 


IV 



3.10 Velocity profiles at the pipe midlength at select times during the flow 

cycle for Va = 100 

3.11 Friction coefficient variation at the pipe midlength during one half 

flow cycle for select Va 

3.12 Streamlines within the left end region at select times during the flow 


cycle for Va = 1 

3.13 Streamlines within the left end region at select times during the flow 

cycle for Va = 30 

3.14 Streamlines within the left end region at select times during the flow 


. 55 


cycle for Va = 100. 


4.1 Velocity profile for fully developed pipe flow plotted in wall coordinates. 78 

4.2 Comparison of experimental and numerical velocity profiles for fully 

developed pipe flow 79 

4.3 Velocity profiles at the pipe midlength at select times during the flow 

cycle for Va = 40 

4.4 Velocity profiles at the pipe midlength at select times during the flow 

cycle for Va — 60 

4.5 Velocity profiles at the pipe midlength at select times during the flow 

cycle for Va = 80 ^ 

4.6 Velocity profiles at the pipe midlength at select times during the flow 

cycle for Va = 100 ^ 

4.7 Velocity profile in wall coordinates at the pipe midlength for Va = 40. 84 

4.8 Velocity profile in wall coordinates at the pipe midlength for Va = 60. 85 

4.9 Velocity profile in wall coordinates at the pipe midlength for Va = 80. 86 


V 



4.10 Velocity profile in wall coordinates at the pipe midlength for Va = 

100 87 

4.11 Friction coefficient over one half flow cycle at the pipe midlength for 

Va = 40 88 

4.12 Friction coefficient over one half flow cycle at the pipe midlength for 

Va = 60 89 

4.13 Friction coefficient over one half flow cycle at the pipe midlength for 

Va = 80 9° 

4.14 Friction coefficient over one half flow cycle at the pipe midlength for 

Va = 100 91 

4.15 Variation of an acceleration source term for the dissipation equation, 

proposed by Koehler (1990) 92 

4.16 Streamlines within the left end region at select times during the flow 

cycle for Va = 40 93 

4.17 Streamlines within the left end region at select times during the flow 

cycle for Va = 60 94 

4.18 Streamlines within the left end region at select times during the flow 

cycle for V a = 80 95 

4.19 Streamlines within the left end region at select times during the flow 

cycle for Va = 100 96 

4.20 Grid for nozzle-ended model 97 

4.21 Friction coefficient during the first half of the flow cycle for the nozzle- 

ended model 98 

4.22 Axial velocity profiles obtained from the Lam-Bremhorst and Kim 

turbulence models 99 


vi 



4.23 Normalized turbulent viscosity profiles fit/p obtained from the Lam- 

Bremhorst and Kim turbulence models 100 

4.24 Normalized turbulent kinetic energy profiles 2k/U^ >max obtained from 

the Lam-Bremhorst and Kim turbulence models 101 

4.25 Normalized dissipation rate profiles eu/U^ max obtained from the 

Lam-Bremhorst and Kim turbulence models 102 

4.26 Axial velocity profiles in wall coordinates for the Lam-Bremhorst 

and Kim turbulence models. Experimental data and universal profile 
shown for reference 103 

4.27 Friction coefficient over one half flow cycle at the pipe midlength for 
the Lam-Bremhorst and Kim turbulence models. Experimental data 

and steady flow results shown for reference 104 

4.28 Turbulent viscosity within smooth nozzle end at flow reversal. ... 105 

4.29 Turbulent viscosity within smooth nozzle end 2.5° after flow reversal. 106 

4.30 Turbulent viscosity within smooth nozzle end 5.7° after flow reversal. 107 

4.31 Turbulent viscosity within smooth nozzle end 9.2° after flow reversal. 108 

4.32 Turbulent viscosity within smooth nozzle end 13.0° after flow reversal. 109 

4.33 Turbulent viscosity within smooth nozzle end 17.1° after flow reversal. 110 

4.34 Normalized turbulent kinetic energy behind turbulence generating 

grid for Kim model m 

4.35 Normalized turbulent kinetic energy behind turbulence generating 

grid for Lam-Bremhorst model 112 

4.36 Variation of the RNG coefficient with t? 113 

4.37 Profile of T] for Re = 12000, based on the Lam-Bremhorst data set. 114 

4.38 Profile of the RNG coefficient for Re = 12000, based on the Lam- 

Bremhorst data set 115 


Vll 



4.39 Comparison of the axial velocity profiles of the RNG and Lam- 

Bremhorst models for Re = 12000 116 

4.40 Comparison of the turbulent viscosity profiles of the RNG and Lam- 

Bremhorst models for Re = 12000 117 

4.41 Comparison of the kinetic energy profiles of the RNG and Lam- 

Bremhorst models for Re = 12000 118 

4.42 Comparison of the dissipation profiles of the RNG and Lam-Bremhorst 

models for Re = 12000 119 

4.43 Comparison of the wall coordinate velocity profiles of the RNG and 

Lam-Bremhorst models for Re = 12000 120 

4.44 Comparison of the axial velocity profiles of the RNG and Lam- 

Bremhorst models for Re = 100000 121 

4.45 Comparison of the turbulent viscosity profiles of the RNG and Lam- 

Bremhorst models for Re = 100000 122 

4.46 Comparison of the kinetic energy profiles of the RNG and Lam- 

Bremhorst models for Re = 100000 123 

4.47 Comparison of the dissipation profiles of the RNG and Lam-Bremhorst 

models for Re = 100000 124 

4.48 Near-wall comparison of the dissipation profiles of the RNG and Lam- 

Bremhorst models for Re = 100000 125 

4.49 Comparison of the wall coordinate velocity profiles of the RNG and 

Lam-Bremhorst models for Re = 100000 126 

4.50 Friction coefficient over one half flow cycle at the pipe midlength for 

RNG model with Re max = 12000, Va = 80 127 


Vlll 



5.1 

5.2 

5.3 

5.4 

5.5 

5.6 

5.7 

5.8 

5.9 

5.10 

5.11 

5.12 

5.13 

5.14 

5.15 


Development of the Nusselt number along the pipe length for a fixed 

wall temperature boundary condition 140 

Development of the Nusselt number along the pipe length for a pre- 
scribed flux boundary condition 141 

Temperature profiles at the axial middle of the pipe for Va = 1. . . 142 

Temperature profiles at the axial middle of the pipe for Va — 30. . 143 

Temperature profiles at the axial middle of the pipe for Va = 100. . 144 

Temperature contours within pipe during first half cycle for Va = 1. 145 

Temperature contours within pipe during first half cycle for Va = 30. 146 
Temperature contours within pipe during first half cycle for Va 100.147 
Cycle- averaged dimensionless bulk temperature for Va — 1, 30 and 

100 148 

Nusselt number versus axial position during the first half cycle for 

Va — 1 149 

Nusselt number versus axial position during the first half cycle for 

Va = 30. 150 

Nusselt number versus axial position during the first half cycle for 

Va = 100 151 

Axial- averaged Nusselt number variation during the cycle for Va 1. 152 

Axial- averaged Nusselt number variation during the cycle for Va 30.153 
Axial- averaged Nusselt number variation during the cycle for Va = 

100 154 


5.16 Phase comparison between heat transfer quantities at the axial mid- 
dle of the pipe for Va = 1 

5.17 Phase comparison between heat transfer quantities at the axial mid- 
dle of the pipe for Va = 30 


IX 



5.18 Phase comparison between heat transfer quantities at the axial mid- 
dle of the pipe for Va = 100 157 

5.19 Phase comparison between heat transfer quantities at x/D = 10.4 for 

Va = 1 158 

5.20 Axial-averaged Nusselt number based on T w - T in for Va=l. ... 159 

5.21 Axial-averaged Nusselt number based on T w - Tj„ for Va — 30. . . . 160 

5.22 Axial-averaged Nusselt number based on T w - Tin for Va — 100. . . 161 

5.23 Amplitude, mean and phase shift for the laminar Nusselt number 

correlation as a function of Va 162 

6.1 Nusselt number in developing turbulent flow 173 

6.2 Temperature profiles at the axial middle of the pipe for Va — 40. 174 

6.3 Temperature profiles at the axial middle of the pipe for Va = 60. 175 

6.4 Temperature profiles at the axial middle of the pipe for Va = 80. . 176 

6.5 Temperature profiles at the axial middle of the pipe for Va = 100. . 177 

6.6 Temperature contours within pipe during first half cycle for Va — 40. 178 

6.7 Temperature contours within pipe during first half cycle for Va = 60. 179 

6.8 Temperature contours within pipe during first half cycle for Va = 80. 180 

6.9 Temperature contours within pipe during first half cycle for Va = 100.181 

6.10 Cycle-averaged dimensionless bulk temperature for Va = 40, 60, 80 

and 100 182 

6.11 Nusselt number versus axial position during the first half cycle for 

Va = 40 183 

6.12 Nusselt number versus axial position during the first half cycle for 

Va = 60 184 


x 



185 


6.13 Nusselt number versus axial position during the first half cycle for 

Va = 80 

6.14 Nusselt number versus axial position during the first half cycle for 


Va = 100 iOD 

6.15 Axial-averaged Nusselt number variation during the cycle for Va = 40.187 

6.16 Axial-averaged Nusselt number variation during the cycle for Va 60.188 

6.17 Axial-averaged Nusselt number variation during the cycle for Va = 80.189 

6.18 Axial-averaged Nusselt number variation during the cycle for Va = 


6.19 Phase comparison between heat transfer quantities at the axial mid- 
dle of the pipe for V a = 40 

6.20 Phase comparison between heat transfer quantities at the axial mid- 
dle of the pipe for Va = 60 


6.21 Phase comparison between heat transfer quantities at the axial mid- 
dle of the pipe for Va = 80 

6.22 Phase comparison between heat transfer quantities at the axial mid- 
dle of the pipe for Va = 100 

6.23 Phase comparison between test functions 

6.24 Turbulent kinetic energy profiles at the axial middle of the pipe dur- 
ing transition, for Va = 80 

6.25 Axial-averaged Nusselt number based on T w - T in for V a = 40. . . . 

6.26 Axial-averaged Nusselt number based on T w - T tn for Va = 60. . . . 

6.27 Axial-averaged Nusselt number based on T w - T in for Va = 80. . . . 

6.28 Axial-averaged Nusselt number based on T w - Ti n for Va = 100. . . 


193 

194 

195 

196 

197 

198 

199 

200 


6.29 Amplitude, mean and phase shift for the turbulent Nusselt number 


correlation as a function of V a 


201 


xi 



NOMENCLATURE 


A r 

amplitude ratio, or 

C f 

skin friction coefficient 

°P 

specific heat at constant pressure 


turbulence model constants 

D 

pipe diameter 

fn, fn h 

low Reynolds number functions 

G 

production term in turbulence transport equations 

h 

heat transfer coefficient 

hp 

horsepower 

i 

V=T 

Jo 

Bessel function of the first kind, order zero 

k 

turbulent kinetic energy or thermal conductivity 

kt 

turbulent thermal conductivity 

£ 

mixing length 

L 

pipe length 

Nu 

Nusselt number 

P 

mean pressure 

Pr 

molecular Prandtl number, ^ 

Pr t 

turbulent Prandtl number, 

q 

heat flux 

r 

radial coordinate direction 

R 

pipe radius, coefficient in RNG k-t model 

Re 

Reynolds number, 

Re t 

turbulent Reynolds number, 


Xll 



p 


density 

a k , <r e Prandtl numbers for k, t 

t shear stress or nondimensional time 

<f> phase shift 

X nondimensional length 

u> angular velocity of the bulk flow oscillation 

Superscripts 

* redefined value 

' fluctuating part 

— average quantity 


Subscripts 


b bulk value 

c critical value for transition 

fd fully developed 

coordinate directions in index notation 
in inlet conditions 

m mean quantity 

max maximum value 

t turbulent quantity 

w value at the wall 

0 limiting value 


xiv 



Re y 

turbulent Reynolds number, 

S 

mean rate of strain 

t 

time 

T 

temperature 

77 

turbulence intensity 

u 

mean axial velocity component 

u+ 

dimensionless velocity, — 

U T 

friction velocity, 

u b 

bulk flow velocity 

V 

mean radial velocity component 

Va 

dimensionless frequency of oscillation, e!i — 

Xi 

axial coordinate direction 

Xmax 

amplitude of bulk fluid axial displacement, — 

y 

distance from the wall 

y + 

dimensionless distance from the wall, ^ 

Greek symbols 


Kronecker delta function 

e 

dissipation rate 

V 

ratio of turbulent to mean strain time scale 

9 

crank angle 

0 

dimensionless temperature 


dynamic viscosity 


turbulent viscosity 

Heff 

effective viscosity, p. + p t 

V 

kinematic viscosity, pj p 


xiii 



Chapter 1 


Introduction 


Unsteady flow is present in man, machine and nature. The flow of blood in arteries 
and capillaries in the human body is pulsatile - composed of a mean flow superposed 
with an oscillating component. The tides that wash in and out of rivers, harbors and 
estuaries are unsteady flows with very long periods of oscillation. Many engineering 
devices employ pulsatile and oscillating flow. Pulsating flow is defined here as a 
periodic flow with a net displacement of fluid over each flow cycle. Oscillating flow 
is defined as a periodic flow with a zero mean over each cycle. The subject of this 
thesis is oscillating flow and heat transfer in pipes which make up the heater and 
cooler sections of the NASA Space Power Research Engine (SPRE) currently under 
development. This engine uses the Stirling cycle as the thermal energy converter 
in a power plant for future space applications. The information presented in this 
thesis will of course be applicable to the design of many types of machinery which 
employ oscillating flow and heat transfer. 

1.1 Invention of the Stirling Engine 

What is known today as the Stirling engine was first proposed by Reverend Robert 
Stirling in his patent application of 1816. This closed air-cycle engine was expected 
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to compete with the steam engine, which was a crude, inefficient and dangerous 
device at that time. Other air-cycle engines had been proposed, but none included 
Stirling’s economizer, referred to nowadays as the regenerator. This single compo- 
nent raised the efficiency of the device to the point where it could run on its own 
power, a feat that many contemporary engine designs could not accomplish. Stirling 
recognized that heat could be captured within the regenerator during expansion, 
and then absorbed back into the working fluid later in the cycle. This physical 
insight is notable, since Carnot’s technique of heat engine analysis would not be 
formalized until 1824. An engine of this design was used for a short time in 1818 to 
drive a quarry pump, before the iron heater vessel failed under thermal stress and 
oxidation. 

After this, Robert Stirling’s attentions were directed increasingly toward the 
church, leaving the younger brother James, himself a trained engineer, to pursue 
further development of the engine. The next significant Stirling engine was built in 
1828, employing a small auxiliary compressor to raise the air pressure and thereby 
increase the efficiency. James noted that the engine produced about 20 hp upon 
startup, but that the power fell off as the engine warmed. Apparently in the absence 
of knowledge of Carnot’s recent work, Stirling reasoned that cooling must be applied 
to the engine in order that the “requisite difference of temperature” be maintained 
on opposite ends of the regenerator. Cooling was not addressed in Robert s original 
patent. A water-cooled engine was built in the 1830’s to test the design. A patent 
for the water-cooling feature was granted in 1840. 

Perhaps the first commercially successful Stirling engine was installed at the 
Dundee Foundry in 1841. The water-cooled engine produced 21 hp with an efficiency 
of about 6.5%, similar to the available steam engines. James had by this time also 
reconfigured the firebox in order to eliminate direct radiant heating of the engine 
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and thereby reduce the thermal stresses imposed on the heater vessel. Uniform 
heating was obtained by carefully directing the flue gases to flow around the heater. 
This engine was unable to satisfy the needs of the growing foundry, and was replaced 
in 1843 by a double-acting version of 45 hp and 18% efficiency. No steam engine 
of the time, nor any air-cycle engine of the nineteenth century could compete with 
this efficiency. Nevertheless, the engine was replaced by a steam engine some three 
years later after a series of heater failures which idled the foundry on three separate 
occasions. 

By the middle of the nineteenth century the steam engine had become a reliable 
and reasonably efficient engine to power the industrial revolution. Most large power 
requirements during this period were met by steam engines. Many low power closed 
air-cycle engines were proposed and offered commercially during the end of the 
nineteenth century and running well into the twentieth. At low power, these engines 
could compete with steam in terms of efficiency (typically on the order of 1%), and 
were much simpler and cheaper to build and operate at a time when electricity 
was not widely available. These small engines may well have kept the spirit of the 
Stirling engine design alive during a period when other engines obtained dominance. 

1.2 Renaissance of the Stirling Engine 

The steady improvement of the steam engine, the development of reliable internal 
combustion engines and the widening distribution of electric power eliminated all 
commercial Stirling-type engines by the mid 1930’s. Shortly thereafter the design 
was rekindled by a company seeking to increase the market for it’s own radio equip- 
ment. Philips Company of the Netherlands sought to devise a power source that 
would enable them to sell radios in parts of the world where electric power was not 
yet available over a conventional distribution grid. The external combustion Stirling 
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engine was considered as a means of providing electric power and was preferred over 
storage batteries, which would eventually require recharging and were heavy and 
difficult to transport. The Stirling engine offered several advantages over other heat 
engines, including quiet operation and the ability to run on poor fuels. The Philips 
team was able to build and run a 16 watt engine by 1938. This was the first in a 
long series of engines built by or in association with Philips. By 1950, Philips was 
able to manufacture 150 generator sets of 200 watt capacity. This was closest any 
of the engines came to commercial success. (On a related note, Stirling refrigerators 
developed by Philips were commercially successful for cryogenic cooling and gas liq- 
uefaction.) The Philips engines were tested in household fans, boats, self-propelled 
lawn mowers, cars and buses. An artificial heart pump and several stationary en- 
gines were also tested. The pistons of the Philips engine designs were connected by 
linkages, being referred to as kinematic engines. Philips abandoned Stirling engine 
research in 1979. Theirs was by fax the longest development effort. A history of 
the Philips Stirling engine projects, along with the early history of Stirling engines, 
is given in the very thorough text of Hargreaves (1991), from which much of the 
preceding discussion has been drawn. 

The Philips kinematic engines were never commercially successful for a number 
of reasons, including 

• inadequate lubrication of the pistons 

• gradual fouling of the regenerator 

• inadequate life of the seals used in various locations 

• poor heat transfer into the hot end of the engine 
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These issues highlight the difficulties involved in bringing the otherwise simple 
Stirling engine into common use. 

1.3 The Stirling Engine in Space 

The space race began in 1957 with the launch of Sputnik, providing further mo- 
tivation to develop the Stirling engine. General Motors and Philips had by 1959 
produced a 3 kW generator set driven by a solar-powered Stirling engine. The 
hot end of the engine would be heated by liquid metal circulating through a solar 
collector array. The single- cylinder engine used a rhombic drive mechanism and a 
centrifugal oil circulation system to provide lubrication under microgravity. The en- 
gine met the design power specification, but was overshadowed by the development 
of photo-voltaic cells, which provided sufficient power to satisfy the newly available 
low-power solid-state electronics. 

The manned space station project being developed by NASA will require sig- 
nificantly more power than the simple satellites of the sixties. The engine used 
for this space power application must not suffer from any of the limitations of the 
kinematic engine, and must have a considerable useful fife with (preferably) zero 
maintenance. One means of eliminating some of the problems mentioned above 
is to use a free-piston design. In this design the displacer and power pistons are 
not rigidly connected, but instead communicate with one another through pressure 
forces transmitted by the working fluid. An advantage of this configuration is the 
tendency for the frequency of operation to be fixed by the natural frequency of the 
displacer piston, a parameter which the engine designer controls. Varying load lev- 
els are satisfied in a very natural way by self-induced changes in piston strokes and 
phase angle (West, 1986). As a free-piston engine, journal bearings are not used, so 
mechanical friction occurs only between the pistons and cylinders. The “loading” 
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on these surfaces is symmetrical and negligible in the microgravity environment. As 
a result, the engine requires no lubrication (other than that provided by the working 
fluid). Fouling problems are thereby also eliminated. The entire volume of working 
fluid can be hermetically sealed into the engine casing such that reciprocating and 
rotating shaft seals are eliminated. Insuring adequate heat transfer from the heat 
source to the working fluid remains as a design challenge, providing motivation for 
the present research. Details of the mathematical models used to design free-piston 
Stirling engines axe found in Walker and Senft (1985). An overview of the Stirling 
SPRE development program is given by Dochat (1990). Recent SPRE test results 
obtained at NASA Lewis Research Center are discussed in Cairelli, et al. (1991). 
Optimization of the SPRE design is taking place in a somewhat empirical fashion, 
with adjustments to various components being considered or tested (Wong, et al. 
1992). Improved understanding of the flow and heat transfer in oscillating flow 
would mi ni miz e these expensive hardware design iterations. 

1.4 Other Stirling Engine Applications 

Several ongoing Stirling engine development efforts will be discussed briefly in this 
section. They give some indication of the breadth of applications of the Stirling 
engine. 

Detroit Diesel Corporation and Stirling Thermal Motors Inc. are jointly build- 
ing a 4-cylinder kinematic engine producing 25 kW at 1800 RPM (Bennethum, et 
al. 1991). The engine uses a variable angle swashplate drive to provide power 
control. The crankcase is pressurized to reduce the demand on the reciprocating 
shaft seals. A more conventional and reliable rotating seal is used on the output 
shaft. The engine proposes to solve many of the interminable problems associated 
with kinematic Stirling engines, including lubrication, regenerator contamination 
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and seal life. Cost control efforts have resulted in an engine which is anticipated 
to be competitive with diesel engines of similar size. The modular design allows 
the same power unit to be applied to portable and marine generator sets, hybrid 
electric vehicles, solar conversion and refrigeration. 

Free piston Stirling engines are also being considered for earth-bound applica- 
tions. Mechanical Technology Inc. has considerable experience with free-piston 
engines (Goldwater, 1990). These engines are being considered for residential co- 
generation, heat pumps, unmanned submersibles and hybrid electric vehicles. In the 
later application a low-powered Stirling engine would drive a sealed linear generator 
and act as a range extender for a battery powered vehicle. The size, weight and 
cost of available generators has fallen significantly in recent years. De Graaff (1991) 
argues for the use of a free-piston Stirling engine and sealed linear generator set as 
the sole power source for electric vehicles, given the weight of currently available 
storage batteries. 

Stirling Technology Company, Sandia National Laboratory and others have co- 
operated to design a 25 kW solar power converter system (Wallace, et al. 1991). 
This system uses a free-piston Stirling engine, but drives a conventional rotary in- 
duction generator. The power pistons of this engine act directly as a hydraulic 
pump. A separate hydraulic motor couples the engine to the generator. Besides the 
Stirling engine and collector, most components of the system are commercially avail- 
able products. The design operating life of the Stirling engine itself is 60000 hours, 
a figure presently possible only through the use of the free-piston configuration. 

The free-piston Stirling engine can also be configured with “pistons” of water 
(Fauvel, et al. 1990). Such an engine is referred to as a Fluidyne. In this arrange- 
ment, the addition of heat to the Stirling cycle can be used to pump water. This 
system can be applied in irrigation systems in areas of the world where conventional 
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pumping engines and fuels are prohibitively expensive. 

Many of the papers referenced above are found in the proceedings of the Interso- 
ciety Energy Conversion Engineering Conference, which meets annually. This forum 
includes Stirling engines and related technologies, and is often the first publication 
to present these papers. 

1.5 Literature Survey of Oscillating Flow and 
Heat Transfer 

This thesis considers methods of improving the modeling of friction and heat trans- 
fer in Stirling engine performance codes. The literature survey will be restricted 
to those papers which offer guidance for modeling flow and heat transfer under 
conditions which are similar to those in the heater and cooler sections of Stirling 
engines. Turbulence modeling has been shown to be a possible limiting feature of 
heat transfer predictions and will receive close attention before heat transfer papers 
themselves are reviewed. This literature survey updates the thorough survey of 
Koehler (1990) and extends the scope to include heat transfer in oscillating flow. 

1.5.1 Laminar Oscillating Flow 

1.5. 1.1 Experiments 

Fully developed l ami nar oscillating flow was first studied experimentally by Richard- 
son and Tyler (1929), who used a crank-driven piston-cylinder to produce flow of 
nearly sinusoidal bulk variation. They observed that the peak in the velocity profile 
moved towards the pipe wall at sufficiently high rates of oscillation. Gaver and 
Grotberg (1986) measured the axial velocities in a tapered channel with oscillating 
flow and summed the cycle-averaged flow. They found one steady streaming cell 
on each side of the symmetry plane for moderate oscillation rates, and two cells 
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on each side at higher rates. The steady streaming is predicted by theory and is 
observed in acoustic studies. 

1. 5.1.2 Analytical and Numerical Solutions 

The analytic solution corresponding to the Richardson and Tyler flow was formu- 
lated by Sexl (1930). Uchida (1956) evaluated Sexl’s solution for pulsatile (non-zero 
mean) flow. The velocity gradient was used to evaluate the work required to main- 
tain the pulsating flow, and the dissipation within the flow. The dissipation is found 
to be higher than that of steady flow at the same mean flow rate. Drake (1965) 
obtained a series solution using separation of variables for oscillating flow in a rect- 
angular channel. The peak in the axial velocity was found to occur near a corner, 
in agreement with the previously quoted work of Richardson and Tyler. Iguchi, et 
al. (1992) solved the laminar flow in the entrance region of a square duct, finding 
good agreement with their experimental results, and proposed an expression for 
approximating the length of the entrance region. 

1.5.2 Transitional and Turbulent Oscillating Flow 

Until recently, relatively few papers on turbulent oscillating flow have been avail- 
able. Unlike the laminar case, analytic studies are not possible and laboratory 
and numerical experiments are difficult and expensive. Many of the experimental 
studies provide only qualitative results describing the flow patterns and turbulent 
structure. Numerical studies often suffer from a lack of data suitable for validation 
of the models. A handful of experiments have been performed recently to support 
development of the SPRE engine. This effort was deemed necessary to improve the 
predictions of the available engine design codes. These whole-engine design codes 
were timed to duplicate test data. It was originally assumed that the poor perfor- 
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mance of the models was due to underprediction of viscous flow losses (Tew, el al 
1990). The design codes were then tuned by artificially augmenting the friction loss 
coefficients. In fact, adjustments to the codes cannot be indiscriminately applied to 
the flow losses alone: thermodynamic losses also result from seal leakage and heat 
losses through the pistons and cylinders. Moreover, the distribution of these losses 
varies from one engine design to another. Fundamental research on the important 
components that make up the Stirling engine is required in order to accurately 
characterize each source of losses. 

1.5. 2.1 Experiments 

Park and Baird (1970) observed the decay of oscillations in a u-tube manometer 
and established criteria for the critical Reynolds number (Rc ) for transition based 
on the oscillation rate. End-effects were observed to produce transition at a lower 
Rc in short tubes, while flows that were laminar over most of the oscillating flow 
cycle tended to have higher Rc . In the later cases, the laminar flow may not be 
able to respond and transition might not occur before the next flow reversal. Chan 
and Baird (1974) measured the overall rate of energy dissipation in an oscillating 
liquid col umn and deduced the wall friction. The wall shear was significantly higher 
than for the corresponding steady flow laminar boundary layer. They hypothesize 
that the flow is intermittently turbulent, but note that the transition from laminar 
to turbulent flow must be more gradual than for steady flow. The critical Reynolds 
number for oscillating flow in straight pipes was measured by Sergeev (1966). Sim- 
ilar results were obtained by Hino, et al. (1976) and Ohmi, et al. (1982). 

Tu and Ram aprian (1983) studied sinusoidally pulsed turbulent flows at a fairly 
high mean Reynolds number of 50000. ^Vhen the pulsations matched the character- 
istic frequency of the turbulence (3.5 Hz), the time mean flow was found to deviate 
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from steady flow. The turbulence assumed a nearly constant value throughout the 
cycle. At lower pulsation rates the turbulence was nearly absent during acceleration. 
The law of the wall was violated in all cases during some portion of the cycle. 

Seume (1988) performed experiments to observe the behavior of oscillating flow 
during transition. The test rig consisted of a long straight pipe connected to a 
piston-cylinder driven by a scotch yoke mechanism. The scotch yoke provides si- 
nusoidal piston motion. The length-to-diameter ratio {L/D) of the pipe was 60, 
matching the heater tubes of SPRE. The test fluid was atmospheric air. Data on 
fluid velocity and rms fluctuations was taken using hot film anemometry over a va- 
riety of operating conditions. This allowed for quantitative analysis of the structure 
of the flow throughout the flow cycle at several locations within the pipe. Two 
transition mechanisms were identified. In the first mechanism, fluid within the pipe 
remains laminar, but turbulent fluid from outside the pipe is drawn in and swept 
downstream. The passage of this turbulent slug is detected by the instrumentation 
as transition. Once the slug passes, the fluid may or may not revert to laminar flow, 
depending on local conditions. In the second transition mechanism, the boundary 
layer grows in the normal fashion. At a sufficient velocity, transition may occur if 
local instabilities become laxge. Like previous researchers, Seume found that lam- 
inar low prevailed during accelerating portions of the cycle. Transition occurred 
very near the peak flow rate, with turbulent conditions lasting through most of the 
decelerating portions of the cycle. 

Friedman (1991) continued the work of Seume, taking detailed measurements at 
the same operating point as the SPRE heater tubes. Single hot-wire anemometry 
was used to collect data on the mean and rms fluctuations of the axial velocity. 
Cross-wire anemometry was used for the mean and rms fluctuations of the radial 
velocity component and one Reynolds stress component — uV. This data was used 
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to test the ability of the present numerical models to predict the damping effects of 
acceleration and the convection of a turbulent slug, as will be discussed in chapter 
4. 

Fishier and Brodkey (1991) performed visualization experiments on oscillating 
flow under operating parameters similar to those in arteries of humans and large 
animals. (This parametric range includes that of the SPRE engine heaters and 
coolers.) The flow was driven through a long straight pipe by a scotch yoke and 
piston. The experiment employed trichloroethylene as the fluid, which was doped 
with magnesium oxide particles with a maximum diameter of about 20 micron. 
Data collection took the form of high-speed film, shot from a camera which was 
mounted on a sliding carrier driven by the same scotch yoke. The particle motions 
recorded thus represent the deviation from laminar, rectilinear flow. The films were 
studied to obtain qualitative information about the state of the flow. In most cases, 
turbulence was observed only during the decelerating parts of the flow cycle. 

Simon, et al. (1992) have proposed a set of empirical transition criteria for 
oscillating pipe flow. The result is cast in a form suitable for use in one dimensional 
engine simulation codes. The criteria determine when transition will occur based 
on crank angle, distance from the inlet, peak flow rate and oscillation rate. The 
criteria axe based on experimental data. Once transition is predicted to occur, a 
steady flow friction coefficient at the instantaneous flow rate is to be applied. 

1.5. 2. 2 Numerical Solutions 

Koehler (1990) modeled oscillating flow in straight pipes and examined the ability 
of the Lam-Bremhorst low Reynolds number k-e model to predict transition and 
subsequent relaminarization. The model was able to predict both transition and 
relaminarization, but the response was too fast. Koehler used derived inlet bound- 
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ary conditions for the turbulence quantities, and felt that this may have produced 
some of the discrepancy. The use of experimentally derived values was expected to 
provide better performance, although this was left for future work. 

Mankbadi and Mobark (1991) used a standard k-t model with wall functions 
to test the ability of the model to account for unsteady effects. They recommend 
a dimensionless oscillation rate as the criterion for applicability of the turbulence 
model for unsteady flows. They note that a low Reynolds number model may be 
more appropriate than a standard model using wall functions. 

Ahn and Ibrahim (1992) computed oscillating, incompressible flow in a pipe un- 
der conditions similar to SPRE. The standard, high-Reynolds number k-t turbulence 
model of Launder and Spalding (1974) was employed for turbulence closure. The 
high-Reynolds number k-t model was found to be inadequate for transitional cases 
(with modest values of Re max ). The behavior improved for cases with lower rates of 
oscillation (quasi-steady behavior). In the case of high Re max , the results followed 
closely the steady correlation at the corresponding instantaneous Reynolds number. 
This may be an artifact of the wall functions boundary condition. 

Ibrahim (1993) is, at the time of this writing, employing a novel means of mod- 
eling transition. The approach taken is to model the flow as laminar until such time 
as a set of transition criteria are satisfied. A k-t turbulence model is then activated. 
Early results show good agreement with the experimental data of Friedman (1991). 
The transition criteria axe based on the work of Simon, et al. (1992). 

1.5. 2. 3 Turbulence Modeling 

The effects of turbulence largely determine the friction and hence the pressure drop 
within the tubes of the Stirling engine heaters and coolers. From a computational 
perspective, the heat transfer problem cannot be solved until the flow is understood. 
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Attention will thus focus on finding an accurate or at least adequate turbulence 
model of oscillating flow before heat transfer results are obtained. 

For complex flows, it is desirable to employ a turbulence model that requires no 
a priori knowledge of the turbulent structure of the flow field. This rules out the 
use of the mixing-length and one equation turbulence models (to be discussed in 
chapter 2). The next class of turbulence models, in order of complexity, is the two 
equation models. The most common model of this class is the k-t model. Launder 
and Spalding (1972) discuss the early development of this type of turbulence model, 
which employs modeled transport equations for the turbulent kinetic energy and the 
dissipation rate. These two values axe used to devise an “effective” viscosity. In 
their so-called standard form, the k-t models are used only for regions where the flow 
is fully turbulent, with wall functions used to provide a link to the wall. To provide 
more generality, the so-called low Reynolds number version of the k-t model was 
devised, first by Jones and Launder (1972). This model used damping functions 
on the turbulent viscosity and dissipation equation sink term. Additional terms 
were also added to the kinetic energy and dissipation equations to facilitate the 
application of boundary conditions and also to avoid singularities in the dissipation 
equation at the wall. The authors expressed the hope that further development 
would provide a way to avoid such terms, which are not physically based. Many 
subsequent low Reynolds number k-t models continued to use damping functions 
to avoid singularities and force the turbulent viscosity to zero at the wall. Addi- 
tional ter ms in the kinetic energy and dissipation equations are sometimes used in 

combination with the damping functions. 

One model which employs no additional terms is the low Reynolds number 
k-t turbulence model of Lam and Bremhorst (1981). Damping functions are applied 
to the turbulent viscosity and to the sink term of the dissipation equation. Another 
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function is applied to the source term to boost the dissipation rate near the wall 
and hence force the kinetic energy to zero. This and other low Reynolds number 

models are reviewed by Patel, et al. (1985). 

While the low Reynolds number k-e models have been able to predict, at least 
qualitatively, the relaminarization of steady flows subjected to favorable pressure 
gradients, the opposite is true for decelerated flows. These flows are subject to sepa- 
ration, which the models typically underpredict. Hanjalic and Launder (1980) trace 
the problem to an insufficient dissipation rate, giving rise to a too high turbulent vis- 
cosity. This provides excessive damping, often preventing separation. They added 
the irrotational component of the strain rate to the source term which appears in 
both the kinetic energy and dissipation equations. However, the term was preceded 
by a larger coefficient in the dissipation equation. Rodi and Scheuerer (1986) used 
this model and compared its performance with that of the Lam-Bremhorst model 
and a one equation model. The modified model performed well when subjected 
to adverse pressure gradient flows. The one equation model, with an empirically 

prescribed length scale, also performed well. 

Thangam and Speziale (1992) have recomputed the backward facing step prob- 
lem in an effort to find the cause of errors due to the use of the k-e model. They 
found that since the separation point is well defined for this problem, wall func- 
tions could be used without ill effect. They were able to reduce the error in the 
reattachment length to 12% by use of a sufficiently fine grid. Additionally, by in- 
corporating an anisotropic eddy viscosity model, the error was reduced to only 3%. 
The anisotropic treatment appears to add considerably to the complexity of the 
model. 

A recent low Reynolds number k-e model by Yang and Shih (1992) uses only 
one damping function, applied to the the turbulent viscosity. An additional term is 
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added to the dissipation equation, as was done in the Jones and Launder model, and 
the time scale in the dissipation equation is redefined to eliminate the singularity 
as the wail is approached. The time scale k/e is augmented by the addition of the 
Kolmogorov time scale yjv/t which is negligible in the fully turbulent region, but 
has physical significance very near the wall, based on similarity arguments. 

Mansour, et al. (1989) computed turbulent channel flows using direct numerical 
simulation. The results were processed and used as a data base against which 
to compare the performance of low Reynolds number k-t models. Budgets of the 
various terms in the kinetic energy equation demonstrate that the pressure diffusion 
is negligible in comparison to other terms at all locations within the boundary 
layer. The profiles of kinetic energy and dissipation rate are found to be very 
dependent on the form of the damping functions used for the turbulent viscosity 
and the dissipation equation. A modified damping function for the dissipation 
equation source term is tested and found to produce better agreement with the 
direct numerical simulation data. 

Alternate forms of two equation models have been proposed using other variables 
from which to derive the length scale. Wilcox (1988) presented a k- ui model where 
the frequency of turbulent eddies u — t/k is given by a modeled transport equation. 
Speziale, et al. (1990) propose a k-r model, where the turbulent time scale r is the 
inverse of u>. They argue that the transport equation for r can be solved subject to 
natural, as opposed to contrived, boundary conditions, and that the limiting form 
of the turbulent time scale at the wall contains only viscous terms. The resulting 
profile for r is shown to be almost linear near the wall, eliminating the numerical 
stiffness found in the k-t models. 

Multiple time scale models have been proposed by Kim and Chen (1989), Chen 
and Singh (1990), Kim (1990, 1991, 1992), Liou and Shih (1993) and others. These 


16 



models propose to divide the spectrum of energy-containing eddies into two or more 
(usually two) ranges. Transport equations for the kinetic energy and energy transfer 
rate (dissipation) within each scale are solved to emulate the cascade of energy from 
large, energy-containing eddies to the small eddies, where dissipation occurs. 

A new type of k-t model based on Renormalization Group Theory is presented by 
Yakhot and Smith (1992). The method uses spectral averaging to close the transport 
equations for turbulent kinetic energy and dissipation. A new feature of the model 
is an additional source term in the dissipation equation which is responsive to rapid 
strains. The models of Kim (1992) and Yakhot and Smith (1992) will be considered 
in detail in chapter 2. 

1.5.3 Heat Transfer 

1. 5.3.1 Experiments 

Liao, et al. (1985) reported a decrease in the heat transfer coefficient in turbulent 
pulsating flow. Their results conflict with those of some previous workers, but the 
increase or decrease seems to be dependent on frequency and amplitude parameters. 
Available data covers wide ranges on the parameter map. 

Dec and Keller (1989) measured the heat transfer from the wall of a pulse com- 
bustor tailpipe. They observed Nusselt numbers as much as 2.5 times higher than 
those corresponding to steady flow at the same bulk flow rate. The degree of en- 
hancement increased linearly with pulsation frequency. Recently, Dec, et al. (1992) 
reviewed the literature and observed the trends in heat transfer enhancement or 
decrement due to pulsations. They found that flows with small amplitude pulsa- 
tions produced relatively little effect on the heat transfer, with most researchers 
reporting slight decreases. Flows with large amplitude pulsations, such that flow 
reversal occurred periodically, produced enhancements by as much as a factor of 5.0 
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over the corresponding steady flow heat transfer. 

Smith, et al. (1992) are currently building a test rig for oscillating flow, pressure, 
and heat transfer measurements. The rig is designed to cover the wide range of 
operating parameters that occur in Stirling engine heaters, coolers and regenerators. 

Tang and Cheng (1993) used their experimental data and multivariate statisti- 
cal analysis to produce correlations for the cycle-averaged Nusselt number for heat 
transfer between a pipe and oscillating air flow. The correlations give the Nusselt 
number as a function of Re max , V a and A r . The data were taken on a test rig 
consisting of a copper pipe with L/D = 54, with sinusoidally varying bulk flow. 
Data reduction was performed by combining thermocouple measurements of inlet 
and outlet bulk temperatures with the known, uniform heat flux and the first law of 
thermodynamics to produce Nu results without direct temperature gradient mea- 
surements. Variations on the statistical methods are discussed, and one of three 
correlations is recommended on the basis of minimum error and the ability to re- 
produce the steady flow correlation in the limit as Va tends to zero. The results 
cannot be extended to other fluids. 

Heat tramsfer in oscillating flow is being studied experimentally by Simon amd 
Qiu (1993) at the time of this writing. They have found the Nusselt number vari- 
ation over the cycle to be quite complicated. The experimental test geometry is 
similar to that of the present numerical study, but subtle differences prevent direct 
comparison of results. Simon and Qiu use an expansion ratio of 3.33. The ratio 
used in the present numerical work is 2.0. Moreover, the thin brass outer walls of 
the expansion/contraction end regions in the experiment axe cooled by a jacket of 
circulating, chilled water. The flush-square shoulders of the end regions are essen- 
tially adiabatic. The numerical model uses adiabatic outer walls on the expansion 
and contraction regions with heated shoulders. As a result, the experimental data 
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show a pronounced front of cool fluid entering and passing down the pipe. This is 
substantially different than the numerical model, in which the end regions act as 
a store of heated fluid which is drawn into the pipe after flow reversal, producing 
gradual variations of temperature. 

1. 5.3.2 Analytic and Numerical Solutions 

Faghri, et al. (1979) performed an analytic study of heat transfer in periodically 
developed pulsating laminar flow without area change, and found that the heat 
transfer was enhanced by the unsteady convection superposed on the steady flow 
component. The degree of enhancement increased with the frequency of the pulsa- 
tion. Siegel (1987) and Zhang and Kurzweg (1991) computed the axially enhanced 
heat transfer in laminar oscillating channel and pipe flow, respectively. They found 
that the axial heat transfer can be much larger than in heat pipes. This occurs at 
relatively high oscillation rates, so the the degree of convective transport is negligi- 
ble. 

Yakhot, et al. (1987) modeled the heat transfer in turbulent pipe flow using 
a expression for the turbulent Prandtl number Pr t which depends on the ratio of 
turbulent to molecular viscosity. In this way Pr t is allowed to take values which 
vary across the boundary layer. The expression is based on the work of Yakhot 
and Orszag (1986). Good agreement was obtained for fluids with widely varying 
Prandtl numbers. 

Heat transfer from gas to cylinder walls is modeled by Lee (1983) and Jeong 
(1991) under conditions of oscillating pressure. These simplified analytic studies 
consider reciprocating machines where the heat transfer is driven by the tempera- 
ture gradient produced by oscillating gas pressure in a piston- cylinder arrangement. 
The use of a complex Nusselt number correlation is recommended to account for the 
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phase lag between the wall heat flux and the bulk fluid temperature. This mecha- 
nism is different from the heat transfer which occurs due to convection within pipes 
subjected to oscillating flow. Nevertheless, the use of the complex Nusselt number 
for heaters and coolers subjected to oscillating flow has been proposed. This thesis 
will demonstrate that the complex Nusselt number correlation is inappropriate for 
oscillating flow. 

Ibrahim, et al. (1991) noted that heat transfer results differed when using the 
incompressible and thermally expandable fluid assumptions. The difference between 
thermally expandable and compressible models was slight. 

Cho and Hyun (1990) modeled the flow and heat transfer in laminar pulsating 
pipe flow. They found that the cycle-averaged friction coefficient differed little from 
the steady case, but the average Nusselt number behavior was more complex. The 
heat transfer was increased over the steady value for intermediate values of the 
oscillation rate, but decreased for lower and higher rates. 

Ibrahim, et al. (1992) modeled the flow and heat transfer in oscillating pipe 
and channel flows, and over blunt bodies. They found the local friction factors and 
Nusselt numbers were often an order of magnitude higher than the corresponding 
steady flow values, in the vicinity of an area change. The variation of the fluid 
temperature was found to lag behind the flow variation. 
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Chapter 2 

The Mathematical 
Representation of the Problem 


This chapter will provide an overview of the equations used to represent the phys- 
ical problem in a well-posed mathematical form. Laminar flow is considered first. 
Turbulent flow will be considered later. Gas flows typically have negligible body 
forces. In the absence of gravity this condition will be satisfied exactly. 


2.1 Laminar Flow 


Incompressible laminar pipe flow and heat transfer is governed by the axisymmetric 
Navier- Stokes and energy equations 
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Equations (2.1) through (2.4) are solved using the SIMPLER algorithm of Patan- 
kar (1980). Koehler (1990) provides further details about the numerical procedure. 
Results are presented in chapters 3 and 5 for the laminar flow and heat transfer 
solutions, respectively. 

2.2 Turbulent Flow 

Numerical modeling of turbulent flow is a research topic in itself. The purpose of this 
thesis is not to develop new turbulence modeling techniques, but rather to employ 
available methods to study the behavior of oscillating flows. Certain limitations 
on the degree of involvement in turbulent flow modeling must be imposed. For 
the most part, currently available models axe used, with consideration given to the 
possibility of making modest modifications to these models only when warranted. 

The frequency of turbulent fluctuations in the Stirling engine components is far 
higher than the bulk flow oscillation frequency. This allows us to take advantage 
of the Reynolds-averaging procedure, wherein the Navier-Stokes equations are av- 
eraged over a period of time which is long relative to the turbulent fluctuation time 
scale, but shorter than the bulk flow oscillation time scale. This technique provides 
a set of unsteady transport equations for the turbulent flow variables. 

The convection and diffusion of heat dominate in this application, so compression 
work, pressure diffusion and viscous heating effects are not included in the energy 
equation. The variables used below will be considered to be the mean part unless 
otherwise noted. 

Incompressible turbulent pipe flow is governed by the axisymmetric Navier- 
Stokes and energy equations 

^-(uj) = 0 (2.5) 
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The eddy viscosity concept is used for turbulence closure. The turbulent stresses 


for incompressible flow axe modeled as 
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The normal turbulent stresses are conveniently grouped with the pressure. The 
redefined pressure is stated as 

(2.9) 


P = p+-pk 

For the temperature, the turbulent stresses are modeled as 
— -r ptCj, dT 
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( 2 . 10 ) 

Pr t dxj 

The turbulent Prandtl number Pr t is assigned the uniform value of 0.9, which is 
a typical value found in figure 12-9 of Kays and Crawford (1980) and produces 
good results for steady turbulent pipe flow. This treatment simplifies the solution 
of turbulent heat transfer and eliminates the need to specify a variation of Pr t . 
The turbulent Prandtl number does in fact rise as the wall is approached, but the 
turbulent viscosity tends to zero at the same time. 

Once turbulence closure is applied, the governing transport equations will be 
solved in their axisymmetric forms shown below in equations (2.41) through (2.44). 
Means of closing the equations will now be discussed. 

2.2.1 Turbulence Modeling 

2.2. 1.1 The Mixing-Length Model 

The eddy viscosity concept requires a means of specifying the field values of the 
turbulent viscosity p t - Numerous techniques have been proposed for this purpose, 
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beginning with the mixing length theory of Prandtl (discussed in White (1974)). 
This technique assumes that turbulent transport is analogous to molecular trans- 
port, and the correlation is fashioned after the kinetic theory of gases. Kinetic 
theory states that the gas viscosity is proportional to the product of the density, 
mean free path and the acoustic velocity: p oc p£v a . Prandtl sought a similar ex- 
pression for the turbulent viscosity that would produce a turbulent shear expression 
with the same form as the molecular shear. In a simple laminar shear flow, we have 



(2-11) 


For the turbulent stress, the desired form is thus 



( 2 . 12 ) 


The turbulent viscosity was then assumed to take the form 


* = T v (2J3) 

where i is the so-called mixing length. The mixing length was determined empiri- 
cally, and was found to depend on the distance from the wall. 

The mixing length model has proven to be quite reliable in certain well-defined 
flow situations such as fully-developed pipe flow. For more general flow geome- 
tries, specifying t is difficult or impossible. Alternate methods of computing p t are 
required, especially when transport and unsteady effects axe present. 


2. 2. 1.2 The One-Equation Models 

The next level of complexity in turbulence modeling is to generate a transport 
equation for one or more turbulence quantities and express p t as a function of 
these quantities. Depending on the number of transport equations employed, these 
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techniques are referred to as one- or two-equations turbulence models. Prandtl (dis- 
cussed in Schlichting (1979)) derived a one-equation model based on the transport 
equation for turbulent kinetic energy. The turbulent kinetic energy is defined as k = 
+ v^ + w^). Bradshaw, et al. (discussed in Schlichting (1979)) derived a one- 
equation model based on the transport equation for turbulent sheax stress, which 
in turn was derived from the turbulent kinetic energy equation. These transport 
equations can be derived by manipulating the Navier-Stokes equations. 

In Prandtl’s model the turbulent viscosity is considered to be proportional to 
the product of a velocity scale and a length scale 

[it OC U scale T scale • (2.14) 

The velocity scale is taken as 

U scale « k* (2.15) 

The length scale in this case is essentially the mixing length. Bradshaw’s model 
contained a more complicated expression for the length scale which was designed 
specifically for the computation of turbulent boundary layer flow. 

These models are capable of accounting for the convection of turbulent kinetic 
energy, but the difficulty in specifying the length scale for general flows remains. 
The length scale can be inferred from experimental data, but this reduces the scope 
of the model to flows of the same type as the data. A model which requires no data 
is desired. 


2. 2. 1.3 The Two-Equation Models 

A model in which the turbulence length scale is also determined from a trans- 
port equation would be free of the limitations of the mixing length and of the 
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one-equation models. Such models were proposed by Kolmogorov, Harlow and 
Nakayama (1967), and Spalding (all discussed in Launder and Spalding (1972)). 
Kolmogorov proposed a transport equation for the frequency of turbulent fluctua- 
tion k* /£ as a means of specifying the length scale. The model could not be tested 
in 1942 for lack of adequate computing resources. Harlow and Nakayama used the 
transport equation for the dissipation rate k t j£ as the second equation. Spalding 
used a transport equation for the square of the frequency of turbulent fluctuation, 
so the method was similar to that of Kolmogorov. Like Harlow and Nakayama, 
most subsequent use of two-equation models has focused on the use of the dissi- 
pation rate to represent the length scale. This is especially convenient since e is 
required in the sink term of the kinetic energy equation. These models are referred 
to as k-t turbulence models. The transport equations for k and t (or any of the 
other variables) can be derived from the Navier-Stokes equations. The procedure is 
detailed in McComb (1990). 

It is worth noting that the dissipation is defined as 
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The early k-t models considered only those regions of the flow where the turbulence 
was isotropic. The isotropic part of the dissipation is 
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This distinction is of no consequence except near walls where damping occurs pref- 
erentially, producing anisotropy. This will have some consequence when attempts 
are made to extend the model to near-wall regions. 

The length scale can be determined once the turbulent kinetic energy dissipation 
rate t is known. The length scale is given by 


Lacale ^ 

t 


(2.18) 
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Just as in Prandtl’s model, the turbulent viscosity is 


The first successful k-t models were capable of computing the turbulence vari- 
ables only in the equilibrium regions, that is, where the generation and destruction 
rates of turbulence are equal. These are also referred to as the fully turbulent re- 
gions. Near-wall regions which are not in equilibrium are the viscous sublayer and 
the Van Driest layer (discussed in Schlichting (1979)). In the viscous sublayer, the 
presence of the wall damps the turbulent fluctuations to the extent that turbulent 
transport is negligible relative to laminar diffusion. In the Van Driest layer, the 
turbulence is actually intermittent, but the time averaged turbulence may be ap- 
proximated by an exponential decay of turbulence as the edge of the viscous sublayer 
is approached. The so-called wall functions were used to bridge the viscous sub- 
layer and the Van Driest layer. This technique is essentially a means of projecting 
boundary conditions out to the fully turbulent region where the transport equations 
are valid. Wall functions are discussed at length in Patankar and Spalding (1967). 

Wall functions are not strictly applicable for low-speed, unsteady or separated 
flows (Patel, et al. 1985). For oscillating flow, the velocity profile varies rapidly 
with time near flow reversal. A turbulence model which depends on wall functions is 
clearly unsuited for this situation. This shortcoming of the k-t model was recognized 
early on, and efforts to extend the model were made by several investigators. The 
general approach was to modify the existing k-t model, now referred to as the High 
Reynolds Number (HRN) or Standard k-t model. The goal was to modify the model 
in such a way that the governing transport equations could be integrated directly to 
the wall, eliminating the need for wall functions entirely. This family of models is 
referred to as the Low Reynolds Number (LRN) models. Several LRN k-e turbulence 
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models axe reviewed in Patel, et al. (1985). 

Previous work on oscillating pipe flow performed by Koehler (1990) showed 
that the LRN k-e model of Lam and Bremhorst (1981) was capable of predicting 
the turbulent structure in oscillating flow. This model was used extensively in the 
present research, so a detailed description will be given. The results obtained with 
the model will be presented later in this work. 

The Lam-Bremhorst model provides low Reynolds number functions (/i, / 2 ) 
that allow the dissipation equation to be integrated to the wall. This eliminates the 
need to specify wall functions to bridge the viscous sublayer. The turbulent viscosity 
is damped (/ M ) so that it approaches zero at solid walls. The modeled turbulent 
kinetic energy and dissipation rate transport equations used in this model are 



+ pG — pt 


( 2 . 20 ) 


| (O' + + c ' hpG l ~ C2f2P T (2 ' 21) 

The production term G is ^5, where the strain rate in axisymmetric coordinates 


The turbulent viscosity is given by 


2 / dn dv 

l dr dx 


( 2 . 22 ) 


k 2 

fit — pfuCn— 


(2.23) 


The following constants are adapted without change from the standard k-e model 
of Launder and Spalding (1974) 


= 0.09 ci = 1.44 c 2 = 1.92 <7* = 1.0 <7^ = 1.3 


(2.24) 
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(2.25) 


The low Reynolds number functions f„, fi and f? are given by 
/ M = (1 — exp(-0.0163R v )) 2 (l + 


/i = 1 + 




(2.26) 


h = 1 “ ex P(~$t) 


(2.27) 


The low Reynolds number functions approach unity in the fully turbulent region, 
so the model is functionally identical to the standard k-t model in that region. Two 
modest modifications, in the form of upper and lower hmits to the /„ function, 
were employed to enhance performance of the model in transitional cases. The 
value of ffi is not allowed to take values above 1.0. This is a physically reasonable 
restriction. Furthermore, when R t is below 200, /„ is not allowed above 0.5 4- 
0.002572*. These modifications are based on the work of Schmidt and Patankar 
(1988), who found that when /„ is allowed to become too large during intermediate 
iterations, turbulence is damped excessively and transition cannot be predicted. 
For this study, a lower limit of 0.02 was also imposed on The selection of 
this value is somewhat arbitrary, but allowing / M to approach zero was found to 
prevent transition. A similar lower limit was imposed by Zhu and Liu (1991), who 
used the Lam-Bremhorst model and found that transition is prevented when / M is 
allowed to become too small. Upon arriving at a converged solution, none of these 
modifications are actually invoked. They simply help prevent the solution from 
running astray at early stages in the iterative procedure. 

The Lam-Bremhorst model has been used successfully in the present work. Some 
deficiencies were noted, as discussed later in this thesis, and it was decided to in- 
vestigate the use of alternate turbulence models. One alternate modeling procedure 


29 



is the multiple time scale approach discussed by Hanjalic, et al. (1980). This pro- 
cedure divides the turbulence energy spectrum into two or more ranges, and solves 
modeled transport equations for the turbulence variables within each range. This 
provides a means of modeling the cascade of turbulence energy from the large scale, 
energy containing eddies to the smallest scales, where viscous dissipation occurs. 
Convenience dictates the use of only two scales. Two-equation turbulence models 
of the k-e type take no account of the cascade effect, and therefore depend more 
heavily on model constants to provide agreement with experimental results. 


2. 2. 1.4 A Multiple Time Scale Model 


A two-scale LRN turbulence model was proposed by Kim (1992), and was applied 
to the present research. In this model, pairs of transport equations are solved for 
each range. The large scale turbulence variables are called the production range 
kinetic energy ( k p ) and the energy transfer rate (tp). The small scale variables are 
the transport range kinetic energy (k t ) and the dissipation rate (e), which represents 
the same dissipation rate as that of the k-e models. The sum of the large and small 
scale kinetic energies is equivalent to the turbulence kinetic energy of the k-e models. 

The k p and ep transport equations used in this model are 



The production term G is the same as that used in the k-e models. 
The k t and e transport equations are 


(2.28) 


(2.29) 
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(2.31) 


d, , , d , . 

(('* + £) J^) + £ ( &,ep2 + w " C,/e<J ) 

The model constants axe given as 

= 0.09 Cj = 0.21 C 2 = 1.32 C 3 = 1.84 c< = 0.32 C 5 = 1.21 
ce = 1.65 a k = 0.75 at = 1.15 ( 2 - 32 ) 

The low Reynolds number functions are given by 

fe p = l -exp(-Ry) (2>33) 


ft = 1 — 0.13exp(— Ry) 


(2.34) 


1 — exp — 02 -Ry — foRy 2 ') 

1 - exp(-fiiRy) 

with 

ft = 0.005 0 2 = 0.001 & = 0.00011 04 = 0.14 


(2.35) 


(2.36) 


The turbulent viscosity is given by 



(2.37) 


where k = ky + k t . 

The performance of this model in oscillating flow was very similar to that of the 
Lam-Bremhorst k-e model. The details are presented in chapter 4. 


2. 2. 1.5 A Renormalization Group Model 

Another turbulence model to consider for possible application to oscillating flow is 
the renormalization group (RNG) k-e model. The derivation of this model begins 
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with the exact transport equations for k and e, but closure is obtained not through 
modeling of the higher order terms, but through the RNG theory, which employs 
spectral averaging to obtain expressions for the turbulent viscosity and the source 
terms of the k and t equations. RNG theory can be considered a model-building 
method for many different types of problems involving lack of exact closure. A 
version of the RNG k-e model which is appropriate for high Reynolds number ap- 
plications is presented by Yakhot and Smith (1992). 

The RNG closure produces three source/sink terms in the e-equation. The 
first two have the same form as those used in the standard k-e models; only the 
coefficients differ. Yakhot and Smith consider the RNG-based closure for these two 
terms to be exact. This is a striking result, and a tribute to earlier workers who 
had only their intuition to guide them in selecting the proper form for closure. 

The new term in the RNG-based e-equation has no formal closure, but Yakhot 
and Smith (1992) propose a form which contains one adjustable constant. This 
term is responsible for rapid-strain response and has no counterpart in the standard 
k-e models. 

The RNG ifc-equation takes the same form as in the standard k-e models. The 
RNG e-equation is: 


l(pe) + A(«€) = A ((, 


dt 




>£) 


e e 2 c m 7 / 3 (* *) 

+cipG- - c 2 p- - 


(2.38) 

1 + 0t) 3 k 

The ratio of turbulent to mean strain time scales tj = Sk/e. The fixed point (lim- 
iting) ratio for high Reynolds number is rj 0 = 4.38. The value (3 as 0.012 is tuned 
to give good agreement with the experimental value of the von Karman constant. 


The other coefficients are 


c„ = 0.0845 ci = 1.42 c 2 = 1.68 <r k = 0.719 <r e = 0.719 (2.39) 
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The values of , Cj, C 2 , <7 * , <7 g and 170 are <dl established by means of the RNG 
theory. Only 0 is adjusted to obtain agreement with experimental results. The 
turbulent viscosity is given by 


k 2 

Ht = pc^— 


(2.40) 


The RNG k-e model was tested for steady pipe flow in the high Reynolds number 
form presented above. The model produced an improvement in the radial distribu- 
tion of turbulent viscosity, relative to the Lam-Bremhorst model, but the difference 
was slight. The model was not applied to oscillating flow due to the limitations of 
the wall function method discussed in section 2.2.1.3 above. 

Modifications to create a low Reynolds number RNG k-e model are discussed 
briefly by Yakhot and Orszag (1986). This author was unable to find a complete 
presentation of the LRN RNG k-e model in the literature. Instead, the HRN version 
discussed above was modified with the addition of the /1 and fi low Reynolds 
number functions from the Lam-Bremhorst model to the first two source terms of the 
e-equation. The last term was modified slightly to prevent unnatural behavior. The 
function fp from Lam-Bremhorst was added to the turbulent viscosity expression. 
The details of these modifications and the results will be discussed in detail in 


chapter 4. 


2. 2. 1.6 Other Turbulence Models 

Turbulence is in fact a convective process. The eddy viscosity models take account 
of the effect of turbulence on the mean flow by attributing the turbulent transport to 
the “turbulent viscosity.” In other words, the convective turbulent transport is mod- 
eled as an enhanced diffusive transport. It is believed that turbulent flow behaves 
according to the Navier-Stokes equations (2.1) through (2.4): the same equations 
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used to model laminar flow. The differences between solutions for laminar and tur- 
bulent flows using these equations are the small spatial scales and unsteady nature 
of turbulent flow. Resolving turbulent flow at the scale of the smallest structures is 
possible, but only at great cost. This type of turbulent flow computation is usually 
referred to as Direct Numerical Simulation (DNS). The required grid spacings and 
time steps used in DNS are so small as to make routine computations of turbu- 
lent flow impractical. The method has been of great utility in producing data sets 
for subsequent statistical analysis, which has led to a better understanding of the 
structure of turbulence. The review paper of Moin (1990) provides a discussion of 
the DNS method and the related Large Eddy Simulation method. 

The Reynolds-averaged Navier-Stokes equations, which are closed by the eddy 
viscosity treatment in this thesis, can also be closed by solving transport equations 
for the second order stress terms of equations (2.6) and (2.7). However, transport 
equations for the second order terms contain third order correlations: this so-called 
turbulence closure problem requires that we truncate, at some level, our effort to 
derive stress transport equations. The Reynolds Stress models (RSM) axe obtained 
by expressing the third order correlations in terms of the mean flow field. This 
also provides a natural means of accounting for the anisotropy of complex turbulent 
flows. RSM models have been successful in capturing some flow phenomena that are 
missed by the eddy viscosity models. Notable examples are the secondary flow in 
square ducts and flows with rotation. However, the number of transport equations 
to be solved is increased significantly relative to eddy viscosity models, even when 
algebraic truncations are applied. There also exists the issue of whether to apply 
wadi functions or integrate the governing equations to the solid walls. Launder 
(1989) points out that the development and testing of LRN stress models lags well 
behind that of the eddy viscosity models. For these reasons, RSM models have not 
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been used for this thesis research. 


2.2.2 The Mean Flow 


Once the field values of fi e jj (= fi + Ht) are available, the time averaged Navier- 
Stokes and energy equations can be solved. The incompressible, axisymmetric forms 
of these equations are 

|H + I|-(rv) = 0 (2-41) 

ox r or 



Equations (2.41) through (2.43) are solved sequentially along with the turbulence 
transport equations using the SIMPLER algorithm of Patankar (1980). Coupling 
enhancement is employed to speed the convergence of the solution of the turbulence 


transport equations. Coupling enhancement is a minor modification to the sequen- 
tial solution procedure which solves the pair of strongly coupled k- and e-transport 
equations several times before proceeding to solve the next variable in the sequence. 
Since constant fluid properties are assumed, the heat equation (2.44) is solved after 


a converged flow solution is obtained. Details of the solution procedure are found 
in Koehler (1990). 
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The oscillating flow and heat transfer problem is solved over the domain shown 
in figure 2.1. Figure 2.1 also shows representative control volumes over which the 
discretized equations were solved (only alternate control volume faces have been 
shown). The grid is refined near the walls of the domain to resolve the steep 
velocity gradients. A somewhat coarser grid is acceptable near the axis and in 
regions remote from the expansion/contraction. 

The inflow boundary conditions are used to create an oscillating flow. The axial 
velocity component at the inflow is specified as 

Ub,in = Ub, ma . x sin(u>t) (2.45) 

This boundary condition is applied to the left end of the domain if sin(ui) is positive, 
and to the right end when it is negative. The radial velocity component is set to 
zero at the inlet. A zero axial gradient boundary condition is used at the outlet end 
for all variables. The no- slip condition is applied to the velocity components at all 
walls. A symmetry condition is used at the axis for all variables. 
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Figure 2.1: Axisymmetric domain geometry and computational grid for pipe with 
expansion/contraction ends. 
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Chapter 3 

The Laminar Flow Solution 


In this chapter the flow solution for laminar pipe flows will be presented. Fully 
developed pipe flow will be defined and some results will be shown to demonstrate 
the concept of fully developed flow and to verify the performance of the computer 
program. Oscillating pipe flow will be considered next. Analytic solutions are 
available un der certain conditions: these will be studied and compared with the 
computer program results. 

3.1 Fully Developed Pipe Flow 

The concept of fully developed flow is used to distinguish between entrance flows 
and flows which are well downstream of a pipe inlet. In engineering flows, the 
entrance region may be only a small fraction of the total pipe length. In these 
cases, piping systems can be designed with reasonable accuracy using pressure drop 
and heat transfer correlations which are based on fully developed flow data. When 
necessary, corrections may be applied later in those cases where entrance effects are 
deemed significant. Incropera and DeWitt (1985) give the length of the entrance 
region (the hydrodynamic development length) for laminar flow as 

f^-) «0.05 Re (3.1) 

\DJ jd 
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Beyond this point the flow is fully developed. The flow is considered to be laminar 
when the Reynolds number is less than about 2300. 

The defining concept of fully developed flow is the absence of axial dependence. 
All terms in the momentum equation which involve axial gradients are thus zero. 
For incompressible laminar flow without body forces, the axial momentum equation 
in axisymmetric coordinates then reduces to 

+*('£) (3 - 2) 
The pressure gradient is a constant by the fully developed flow definition. Equation 

(3.2) can then be integrated twice to obtain the axial velocity profile. By applying 
the no-slip boundary condition at the pipe wall and the zero-gradient condition at 
the axis, the profile may be written as 

„(r) = (if - r 2 ) (3-3) 

The velocity profile is parabolic in the radius. The elliptic flow solver employed 
for the present research must be able to reproduce this result at positions well 
downstream of the pipe inlet. 

Figure 3.1 shows the numerical solution for the fully developed velocity profile 
using three different radial grids. The results were taken at the outlet of a straight 
pipe with a length of 500 diameters. A good solution can be obtained with a 
surprisingly coarse grid. Radial grids with 5, 10 and 20 uniform control volumes were 
used to produce the profiles in figure 3.1. The solution with only 5 control volumes 
is very reasonable, while the solution with 20 control volumes is indistinguishable 
from the exact solution to the resolution of this figure, serving to verify the computer 
program. 
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3.2 Oscillating Pipe Flow 


The study will now be extended to the unsteady case where the bulk flow varies 
sinusoidally. The oscillating flow can be said to assume a periodically developed 
state after a sufficient number of flow cycles have been completed. Fully developed 
oscillating flow was first studied experimentally by Richardson and Tyler (1929). 
They observed that the peak in the velocity profile moved towards the pipe wall at 
sufficiently high rates of oscillation. The governing equations for this flow can be 
solved analytically for positions fax removed from the pipe inlet. This solution is 
discussed below. The corresponding numerical results will then be compared to the 
analytic solution for verification. 


3.2.1 Analytic Solution 


The analytic flow solution corresponding to the Richardson and Tyler study was 
formulated by Sexl (1930). The separation of variables technique was used to obtain 
the axial velocity time- and radial-dependence given by the real part of 


u(r,t) 


K exp iut 
iuj 


1 


J 0 ryJ^iuj/ v 

J 0 R,J^uJv 


(3.4) 


where 


Kexpiut = — — (3.5) 
p dx 

Figures 3.2 through 3.4 show the velocity profiles throughout the oscillating flow 
cycle given by equation (3.4) for Va = 1, 30 and 100. Notice that the Sexl solution 
is pressure gradient-driven. The crankangles of figures 3.2, 3.3 and 3.4 refer to 
the variation of the pressure gradient throughout the cycle. The most significant 
feature of oscillating flow is the increase in the wall velocity gradient, particularly 
at low flow rates. This creates an axial pressure gradient which is higher than that 
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of steady flow at the same Re. The pressure gradient is also out of phase with the 
bulk flow rate. The bulk flow lags the pressure gradient by 90° a sVa becomes large. 
Figure 3.5 shows this phase shift as a function of Va. 

3.2.2 Numerical Solution 

A numerical model used to compute fully developed pipe flow can be extended to 
model unsteady flow by making the boundary flow rate vary with time. For this 
section, the computer program was set up to model a pipe with sudden expansion 
and contraction end regions. The domain geometry and grid are shown in figure 3.6. 
Figure 3.7 compares the axial velocity profiles of the Sexl analysis and the present 
study at the peak instantaneous flow rate with V a — 30, for which the velocity lags 
the pressure gradient by 73°. The agreement between the analytic and finite volume 
solutions is good, particularly near the wall, where the shear stress is computed. 
The numerical solution was obtained with only 16 control volumes across the pipe 
radius, and the profile used in figure 3.7 is taken at the axial middle of the pipe 
section, 25 diameters from either end. 

Figures 3.8 through 3.10 show the axial velocity profile at the middle of the pipe 
at select crank angles during the flow cycle for Va = 1, 30 and 100. These crank 
angles refer to the variation of the bulk flow rate. The velocity profiles for low 
values of Va differ only slightly from the steady flow case, as shown in figure 3.8. 
The profiles axe nearly identical for accelerating and decelerating flow. Compare 
the profiles of figure 3.8 for 58.7 and 121.3° crank angle, for instance. They are 
virtually indistinguishable. The velocity profiles become more complicated at higher 
oscillation rates, as seen in figures 3.9 and 3.10. The effect of oscillation is most 
apparent near flow reversal (crankangles of 180° and 360°). The momentum of the 
core fluid is relatively higher than that of the near-wall fluid, causing a recirculation 
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of fluid within the pipe, even though the bulk flow is instantaneously zero. After 
flow reversal the low-momentum fluid near the wall responds quickly to the axial 
pressure gradient, and accelerates ahead of the core fluid. In figure 3.9 the axial fluid 
does not take the the lead again until just after 90°. The peak in the velocity profile 
is off the axis for most of the flow cycle when the rate of oscillation is high, as seen 
in figure 3.10. In all oscillating cases, as the fluid decelerates the high-momentum 
core fluid flows at a sufficient rate to cause the near-wall fluid to reverse prior to 
bulk flow reversal. This is required by the incompressible continuity equation to 
ensure that the mass rate is constant at all sections along the pipe. At the middle of 
the pipe the next half-cycle is identical to the last, with all flow directions reversed. 

The wall shear is phase shifted in time with respect to the bulk flow velocity 
in oscillating flow. Figure 3.11 shows the friction coefficient as it varies through 
the first half cycle at the midsection of the pipe for Va = 1, 30 and 100. Re max 
is 2000 in each case. The steady flow friction coefficient at the corresponding in- 
stantaneous flow rate is shown for reference. At low rates of oscillation the friction 
coefficient differs little from the steady case (except near flow reversal). As the rate 
of oscillation increases the friction coefficient becomes generally higher during the 
accelerating portion of the cycle and lower during deceleration. During acceleration, 
the low-momentum fluid near the wall responds quickly to the pressure gradient, 
creating an enhanced wall velocity gradient and friction coefficient. During deceler- 
ation, the opposite effect occurs. The near-wall velocity reverses prior to bulk flow 
reversal, and the friction coefficient drops sharply to zero. Figure 3.11 has been pre- 
pared by using the signed values of the near-wall velocity to evaluate the wall shear 
stress. Though unusual, this treatment has been used to illustrate the near-wall 
flow reversal which occurs prior to the bulk flow reversal at 180° crank angle. The 
friction coefficient curves for the oscillating cases all pass through zero when the 
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near-wall flow reverses. In order to obtain the more conventional, unsigned friction 
coefficient, the negative portions of the curves may be simply mirrored across the 
abscissa. 

The flow field is far more complicated within the end regions. Figures 3.12 
through 3.14 show streamline plots in the left end region at select times during the 
cycle for the cases Re max = 2000 and Va = 1, 30 and 100. Flow is to the right for 
crankangles between zero and 180° and to the left for the remainder of the cycle. For 
low values of Va the effects of oscillation are apparent during inflow only near flow 
reversal (see figure 3.12 for Va — 1, 180°), when the momentum of the core fluid 
causes a slight recirculation. In the second half of the flow cycle the momentum of 
the flow is sufficient to cause separation, which in turn traps a large portion of fluid 
within a recirculation zone. This recirculation zone is initiated during outflow and 
persists through a portion of the flow cycle, depending on V"a. As Va is increased, 
less time is available for the growth of the recirculation zone. When the rate of 
oscillation is increased, the recirculation persists over a larger portion of the flow 
cycle. Figure 3.13 shows that the recirculation left over from the previous cycle 
is still relatively strong within the left end region at about 30°. This early in the 
cycle, the incoming flow is forced to move around the recirculating slug of fluid left 
over from the previous cycle. Eventually this slug is sheared apart and by about 
120° the effect of oscillation is barely noticeable. During outflow, the center of the 
separation bubble remains slightly closer to the shoulder of the expansion (where it 
originates) relative to the Va = 1 case. When Va is increased to 100, the effects of 
oscillation are apparent throughout the flow cycle as shown in figure 3.14. During 
outflow the separation bubble is closer yet to the shoulder of the expansion. 
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Figure 3.2: Velocity profiles throughout the oscillating flow cycle for V a = 1 as 
given by Sexl’s analytic solution. 
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Figure 3.4: Velocity profiles throughout the oscillating flow cycle for Va — 100 as 
given by SexTs analytic solution. 
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Phase Shift 



Figure 3.5: Variation of the bulk velocity-pressure gradient phase shift with V a. 
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Figure 3.6: Axisymmetric domain geometry and computational grid for pipe with 
expansion/contraction ends. 
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Figure 3.8: Velocity profiles at the pipe midlength at select times during the flow 
cycle for V a = 1 . 
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Figure 3.9: Velocity profiles at the pipe midlength at select times during the flow 
cycle for V a = 30. 
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Figure 3.11: Friction coefficient variation at the pipe midlength during one half flow 
cycle for select V a. 
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Figure 3.12: Streamlines within the left end region at select times during the flow 
cycle for V a — 1. 
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Crank Angle = 180.0 Crank Angle = 360.0 

Figure 3.13: Streamlines within the left end region at select times during the flow 
cycle for V a = 30. 
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Figure 3.14: Streamlines within the left end region at select times during the flow 
cycle for V a = 100. 
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Chapter 4 

The Turbulent Flow Solution 


In this chapter the flow solution for turbulent flows will be presented. Fully devel- 
oped results for pipe flow will be shown first, as a means of verifying the performance 
of the turbulence models. Results for oscillating flow will then be presented. Com- 
parisons with the available experimental data will also be made, with points of 
agreement and discrepancies noted. 

The operating conditions considered in this and following chapters are, unless 
otherwise noted, representative of the SPRE operating conditions. For SPRE, the 
maximum Reynolds number (Re max ) during the flow cycle is 12000. The oscillation 
rate is given by the Valensi number (Va), which is 80 for the SPRE design point. 

Three different turbulence models were tested during the course of this work. 
Most of the results presented axe based on the Lam-Bremhorst Low Reynolds Num- 
ber (LRN) model. The two time scale model of Kim and the RNG model of Yakhot 
and Smith were also tested. The mathematical bases of these turbulence models 
were discussed in chapter 2. 
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4.1 Fully Developed Pipe Flow 

The definition of fully developed flow does not depend on whether the flow is lam- 
inar or turbulent. The absence of axial variation in all the dependent variables is 
sufficient to satisfy the fully developed condition. The nature of turbulent flow, 
with enhan ced cross stream transport, is such that fully developed conditions are 
typically established within only a few diameters of the pipe inlet. Incropera and 
DeWitt (1985) consider turbulent flow to be fully developed after 10 to 60 pipe 
diameters. 

The increased cross stream transport which occurs in turbulent flow has a dra- 
matic effect on the axial velocity profile. Relative to laminar flow, the velocity 
profile is flatter neax the axis, and the velocity gradient at the wall is much steeper. 
As a result, the pressure gradient required to drive a turbulent flow is much larger 
than for laminar flow. Turbulent flows typically occur in piping systems when the 
Reynolds number is greater than about 2300. 

The velocity profile for turbulent flow has been the subject of a great deal of 
experimental research. The measured profile can be fitted well using the so-called 
law of the wall. The flow can be divided into three regions, depending on the 
distance from the wall. Next to the wall, turbulent fluctuations are damped and the 
molecular diffusion dominates. This is referred to as the viscous sublayer. Far away 
from the wall, the effects of molecular diffusion are negligible in comparison with 
turbulent transport processes. This is called the fully turbulent region. Between 
these layers is a region where the flow is intermittently laminar or turbulent. Time 
averaging reveals that a gradual increase in turbulent activity occurs as this region 
is traversed from the viscous sublayer to the fully turbulent region. This is referred 
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to as the buffer region. For the viscous sublayer, the velocity profile is given by 

u + = y + (4.1) 

where 

u + = — and y + = (4-2) 

u T v 

are the wall coordinates and 

(4.3) 

is the friction velocity. In the fully turbulent region, the profile is fitted to the 
expression 

u + = A log y + 4- B (4.4) 

The form of this expression also suggests the name log-linear for this region. Schetz 
(1984) recommends A = 5.6 and B = 4.9 as the best fit through the available 
experimental data in the log-linear region. Equations (4.1) and (4.4) make up the 
law of the wall. 

The performance of the computer model can be tested against the law of the 
wall profile. Figure 4.1 compares the computed velocity profile (in symbols) and 
the law of the wall (dashed lines). The agreement is good. The computed results 
are based on fully developed turbulent flow using the Lam-Bremhorst turbulence 
model at Re = 50,000. The grid is composed of 42 nominiform grid points across 
the pipe radius, with near-wall grid points clustered tightly together. As is typical, 
the pipe flow profile shows a slight “wake” at the axis. The low Reynolds number 
functions of the Lam-Bremhorst turbulence model do a good job of fairing in the 
viscous sublayer and the log-linear regions. 

Figure 4.2 compares the same predicted fully developed axial velocity profile 
with the experimental data of Laufer (1953). Alternate numerical data points are 
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plotted for clarity. The numerical curve is normalized against it’s own peak value. 
The experimental data is normalized to match the bulk velocity of the numerical 
profile. The agreement is fair. The numerical model produces a somewhat flatter 
profile near the axis, and thus overpredicts the velocity closer to (but not at) the 
wall. This is a common deficiency of the k-e model, due to a slightly overpredicted 
turbulent viscosity near the axis. The numerical model nevertheless produces the 
correct wall shear. 


4.2 Oscillating Pipe Flow 


As was done in chapter 3, the numerical model can be extended to oscillating flow by 
imposing a time-dependent velocity boundary condition. Oscillating flow in straight 
pipes was studied numerically by Koehler (1990). This thesis presents the results 
obtained when the model is extended to the case where expansion and contraction 


regions are applied to the ends of the pipe. 

The motivation for modeling oscillating pipe flow with enlarged end regions is to 
observe what effect the end regions have on the flow within the pipe itself. Koehler 
specified a uniform axial inlet velocity for his numerical model. The boundary 
conditions for the turbulence quantities at the inlet to the straight pipe were derived. 
These boundary conditions were based on the assumption of isotropic turbulence, 
for which 


k m = I ( TI 


(4.5) 


is the definition of the turbulence intensity TI. The available data for fully devel- 
oped turbulent pipe flow was then manipulated to produce an expression for the 
dissipation rate 


— 


l 

TI Re i 



v 


(4.6) 
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This functional dependence on Re follows from the data, but the use of TI is 
somewhat arbitrary. This treatment requires that the flow be allowed several pipe 
diameters in which to develop before the results can be considered representative of 
pipe flow. The same turbulence boundary conditions are used in the present study, 
but the end regions provide a space in which the turbulence can adjust prior to 
reaching the pipe itself. Interesting flow features are also allowed to develop in the 
end regions. For instance, flow separation occurs when high momentum fluid flows 
into an expansion. Upon flow reversal, these features are expected to be swept back 
into the pipe, and their effect on the pipe flow can be observed. 

The results for turbulent flow depend strongly on the values of Re and Va. For 
moderate values of these parameters, the flow is not turbulent throughout the entire 
cycle. Especially for low values of Va, the flow may be laminar for a considerable 
time near flow reversal, when the bulk flow rate is low. As Va increases, turbulent 
effects persist over a larger portion of the flow cycle. In this and following chapters 
we will refer to “turbulent oscillating flow” while recognizing that laminar and 
transitional flow will be present during part of the flow cycle. 

4.2.1 Results for the Sudden Expansion/Contraction Ends 

In this section the results for the pipe fitted with sudden expansion/contraction 
ends are considered. The Lam-Bremhorst turbulence model is used. The results 
are compared with the available experimental data. 

Figures 4.3 through 4.6 show the axial velocity profiles at the pipe midlength 
at select times during the flow cycle for Va = 40, 60, 80 and 100. Re max in each 
case is 12000. For moderate values of Va, the effect of oscillation is apparent only 
near bulk flow reversal. In figure 4.3 the circulation of flow at bulk flow reversal 
is small, but not negligible. Notice, however, that the profiles for crank angles of 
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58.7° and 121.3° axe nearly indistinguishable. When Va is raised, the profiles for 
accelerating flow begin to differ from those of decelerating flow, particularly near 
the wall. Figure 4.6 shows that the profiles for 58.7° and 121.3° are different. The 
profile for 58.7° is somewhat laminar-like. The wall gradient is lower than that of 
the 121.3° profile. For Va = 100, the wall-near fluid at 58.7° has not yet made the 
full transition to turbulence. The profile at 121.3° shows a steeper wall gradient 
since the flow is turbulent at this point in the cycle. This effect of transition is not 
seen in the laminar case, and causes the turbulent flow to be more complicated than 
the laminar flow. 

As was mentioned in chapter 2, the law of the wall is not generally applicable 
for oscillating flow. However, it is nonetheless instructive to observe the predicted 
axial velocity profile in wall coordinates. Figures 4.7 through 4.10 show the velocity 
profile in wall coordinates for Va = 40, 60, 80 and 100. The profiles are measured 
at 90° crank angle at the axial midlength of the pipe. The rate of change of the 
pressure gradient is modest at 90° cra nk angle and the end effects are insignificant 
30 pipe diameters downstream of the inlet. Under these conditions, the velocity 
profile is expected to be at least qualitatively similar to the universal profile, and 
this is born out by the figures. Figure 4.9 also shows the experimental data of 
Friedman (1991) for the SPRE operating condition. These data differ from the 
universal profile, but the differences can be explained as an unsteady effect. Once 
transition has occurred, the velocity profile near the wadi is expected to be steeper 
than for steady flow and the flow near the pipe axis is slightly slower. Figure 4.9 
shows that the experimental velocities are indeed slower near the axis. The Lam- 
Bremhorst turbulence model apparently causes the velocity to follow the universal 
profile too closely, and thereby understates the unsteady effect. 

The variation of the friction coefficient over the first half of the flow cycle is 
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shown in figures 4.11 through 4.14. These results axe taken at the pipe midlength. 
The friction coefficient is defined as 


Cf = 



(4.7) 


This definition of c/ is undefined as the bulk flow rate approaches zero, in which 
case the definition would have no utility. 

The Blasius friction coefficient correlation for steady flow at the corresponding 
instantaneous Reynolds number is also plotted on each figure for reference. Schlicht- 
ing (1979) gives this correlation a 5 


c f = O.miRe- 025 


(4.8) 


The modest Reynolds number dependence produces a rather flat curve over the 
cycle. 

Figure 4.11 shows that the friction coefficient starts high and drops quickly 
during the early part of the flow cycle. This is indicative of the laminar nature 
of the flow at this point in the cycle. At about 25° crank angle the predicted 
friction coefficient drops below the corresponding steady flow value. The flow is 
still l amin ar-like at this point. At 35° the friction coefficient starts to rise up and 
meets the steady flow curve. This behavior is expected at low values of V a , for 
which the flow can be characterized as quasi-steady. As V a increases, the laminar 
behavior of the flow persists farther into the cycle. Figures 4.12, 4.13 and 4.14 
show that transition is predicted at roughly 40°, 45° and 60°, respectively. In each 
case, the flow near the wall reverses prior to bulk flow reversal, causing the friction 
coefficient to drop sharply to zero before rising again and taking a laminar-like 
value. 

The experimental friction coefficient for the SPRE operating condition based 
on Friedman’s data is also shown on figure 4.13. Up to about 40° the agreement 
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between the predicted and measured friction coefficients is excellent. However, the 
experimental data continue to show laminar- like behavior up to nearly 90°. The 
experimental data were measured with a smooth nozzle entrance to the pipe, while 
the numerical results are based on a sudden contraction pipe entrance. Friedman 
(1991) states that the observed transition at this axial location coincides with the 
arrival of a slug of turbulent flow which was generated within the nozzle at the end 
of the previous flow cycle. This slug is then convected back into the pipe during 
the next cycle. This turbulent fluid may be responsible for the fact that the data 
for crank angles of 90° through 130° are above both the predicted and steady flow 
results. The traveling slug issue was investigated further by building a nozzle-ended 
numerical model and subjecting it to the SPRE operating condition. This model 
and the corresponding results will be discussed below. 

A turbulent slug generated within the end region during outflow would not be 
expected to survive the intense shearing action which occurs at a sudden contraction 
after flow reversal. Apart from the traveling slug issue, the experimental data of 
figure 4.13 show that the turbulence within the pipe is damped by the acceleration 
of the fluid up to at least 80° crank angle. The numerical model predicts transition 
at about 45° crank angle, at which time the bulk flow Re is about 8500. The bulk 
flow Re passes 2300 at about 11° crank angle, showing that the model does reflect 
at least some degree of damping due to acceleration. It would be desirable to find a 
means of modifying the model to account for the experimentally observed degree of 
damping. Koehler (1990) proposed a modification to the dissipation rate transport 
equation (2.21) to account for the effects of acceleration and deceleration. The 
modification consists of adding an additional production term which is scaled by 
the instantaneous value of acceleration. The proposed term would take the form 

eK a ( 4 - 9 ) 
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where 


K a 


" -im 


pUi, m „ at 


(4.10) 


and U h is given by U b , max sin(ut). The dissipation rate transport equation would 
then become 


I + 

((M + ft) £) + ^ (O + «**•) - e*4 (4.1D 

The constant c 3 would be fitted against experimental data. The new term would 
thus augment the existing production term G. 

The sinusoidal variation of the bulk velocity in the present research will produce 
an acceleration source term which is proportional to cos(ut). Curve A of figure 4.15 
shows the variation of cos(wt) over the flow cycle. At the beginning of the cycle 
the flow is accelerating. The new source term would provide a positive contribution 
to the dissipation, resulting in the desired damping of turbulence. During the 
second quarter-cycle the contribution would change signs, amplifying turbulence 
during deceleration. However, the new source term contribution would be physically 
unrealistic during the second half cycle. The fluid itself is unaware of the coordinate 
direction. Regardless of the direction of flow, the fluid will sense “acceleration” 
whenever the magnitude of the axial velocity is increasing. In other words, the 
new dissipation rate source term should be positive in the third quarter-cycle and 
negative in the fourth. Curve B of figure 4.15 has the correct sign for the last half 
cycle. 

The proposed modification to the dissipation rate equation given by equation 
(4.10) is physically unrealistic for oscillating flow. The combination of the behavior 
of curve A of figure 4.15 during the first half cycle and curve B during the second 
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would provide the correct sign. However, this treatment results in stepwise dis- 
continuities in the source term contribution at each flow reversal. This is in itself 
physically unrealistic and would require further (ad hoc) modifications to correct. 
The additional modifications would create yet another source of uncertainty in the 
results. It seems that there is no straightforward means of applying an acceleration 
modification to the dissipation rate equation in the case of oscillating flow. This 
treatment will no longer be considered for the present research. 

The streamlines within the left end region are shown in figures 4.16 through 
4.19 for Va = 40, 60, 80 and 100. These figures show that pronounced unsteady 
effects axe present near flow reversal in each case, and the unsteady effect persists 
over a larger fraction of the flow cycle as Va increases. The recirculation bubble 
within the end region during outflow is initiated shortly after 180°. The subplot for 
210° crank angle in figures 4.16 through 4.19 clearly shows the effect of Va. As Va 
increases, less time is available for growth of the recirculation bubble. The bubble 
continues to grow as the flow cycle continues. 

4.2.2 Results for Smooth Nozzle Ends 

In this section we consider the results obtained from the smooth nozzle-ended model. 
The objective of this numerical test is to determine whether the Lam-Bremhorst 
turbulence model can predict the passage of a turbulent slug of fluid entering the 
pipe from the nozzle region. The smooth nozzle-ended pipe is not representative of 
the SPRE design, but was employed during the early stages of experimental oscil- 
lating flow research at the University of Minnesota. Seume (1988) used the smooth 
nozzles in an attempt to prevent large-scale separation at the inlet to the pipe and 
thereby concentrate attention on the fully developed features of oscillating flow. 
This approach was successful except for a small portion of the cycle immediately 
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following flow reversal. The divergence angle of the nozzle was too large to prevent 
separation during outflow. A recirculating bubble of fluid was trapped within the 
nozzle. Upon reversal, this relatively turbulent fluid was convected back into the 
pipe. The fluid entering the test section behind the slug had a low level of turbu- 
lence, since it had passed through a heat exchanger core which anted like a fine-scale 
flow straightener. This heat exchanger was intended to be used during later heat 
transfer experiments, and was included in the flow tests in order to ensure that 
consistent end conditions were employed during both the flow and heat transfer 
phases of the experiments. As the turbulent slug was swept through the pipe the 
instrumentation sensed the arrival at each probe location as transition. Depending 
on the probe location, the fluid would in some cases remain turbulent after the slug 
had. passed, or would revert to laminar flow, depending on local conditions. 

The geometry of the SPRE design will employ different end conditions. One 
end of the heater (and cooler) will be connected to a large plenum, with the other 
end close-coupled to the regenerator matrix. The plenum ends will behave as sud- 
den expansion / contractions . Fluid entering from the plenum will experience intense 
acceleration and shearing. This will break up any large scale flow structures, ren- 
dering the flow in the pipe relatively insensitive to the flow structure in the plenum 
region. Fluid entering the pipes from the regenerator will have a nearly uniform 
axial velocity component, but the turbulence structure will be difficult to specify. 

The traveling slug problem, though not directly applicable to the SPRE design, 
provides an opportunity to test the turbulence model employed for closure of the 
Navier-Stokes equations. Figure 4.20 shows the grid used to model the left nozzle 
end region. The length-to-diameter ratio of the pipe is 60, matching the Seume 
experiment. Each nozzle is 4 pipe diameters in length. The nozzle inlet diameter is 
3.33 times that of the pipe. The nozzle contour is formed by two cubic functions. 
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The nozzle inlet boundary flow is used to create the flow oscillation. 

Figure 4.21 shows the friction coefficient over the first half of the flow cycle, 
plotted along with the experimental data which was also plotted in figure 4.13, at 
the axial middle of the pipe. The results for the nozzle-ended model are essentially 
identical to those obtained with the expansion/contraction model. There is no 
evidence that a turbulent slug is convected into the pipe after flow reversal, or if it 
were it has dissipated before reaching the axial middle of the pipe. 

4.2.2. 1 Shortcomings of the Lam-Bremhorst Model 

The results presented above point out two shortcomings of the Lam-Bremhorst tur- 
bulence model when applied to oscillating flow. First, the model fails to account for 
the full degree of damping which occurs as the flow is accelerated. Figure 4.13 com- 
pared the predicted and measured friction coefficients and showed that the model 
predicted transition too early. For the SPRE operating condition, the flow in the 
pipe passes through Rt — 2300 at about 11° crank angle. The Lam-Bremhorst tur- 
bulence model was able to delay transition to about 45°, but this is still premature 
relative to the experimental results, for which transition did not occur until after 
80° crank angle. The second shortcoming is the inability of the model to predict 
the passage of a turbulent slug of fluid. This feature was demonstrated by the 

experimental results shown in figure 4.21. 

At this stage in the research the two time scale model of Kim (1992) was con- 
sidered to determine whether it would improve the predictions of the acceleration 
damping effect or the traveling turbulent slug. The details of this model where 
discussed in chapter 2. 
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4.2.3 Results of the Kim Model 


In this section the results for the Kim model are compared with those of the Lam- 
Bremhorst model and the experimental data of Friedman (1991). Some general 
comparisons of the predictions of flow and turbulence quantities are first made. 
The transition to turbulence and the acceleration damping behavior of the models 
are compared next. Lastly, the ability of the model to predict a traveling slug of 
fluid ingested from the nozzle end after flow reversal is studied. 

Figures 4.22 through 4.26 compare the results of the Lam-Bremhorst and Kim 
models at peak flow for the Re-max — 12000, Va = 80 case. The sudden expan- 
sion/contraction model results are used, with data taken at the axial middle of the 
pipe. Figure 4.22 shows the axial velocity profiles. The Kim result is not as flat 
as that of the Lam-Bremhorst model, but the wall gradients appear to be similar. 
The Kim velocity profile is slightly better than the Lam-Bremhorst result in that 
the velocity profile is not excessively flat (compare with figure 4.2). 

Figure 4.23 shows the normalized turbulent viscosity profiles for the 
Lam-Bremhorst and Kim models. The Kim result is significantly higher except 
near the wall, where both models show a decay to zero. The Kim result is worse 
than Lam-Bremhorst. According to Schlichting (1979), the turbulent viscosity is 
expected to fall off near the axis, where turbulence production is low. Neither model 
produces a good distribution of turbulent viscosity near the axis. 

Figure 4.24 shows the normalized turbulent kinetic energy profiles for the same 
condition. The Kim model predicts a substantially higher level at all locations 
except the wall. While the boundary condition forces both models to zero at the 
wall, the Kim model produces a far steeper gradient near the wall and a much larger 
overshoot. 
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Figure 4.25 shows the normalized dissipation rate profiles. In this case the Kim 
model produces significantly higher values at all radial positions. 

Figure 4.26 shows the axial velocity profiles of the Lam-Bremhorst and Kim 
models plotted in wall coordinates. The experimental data of Friedman is also 
shown, along with the universal profile in dashed line. The three data sets represent 
the velocity profile at the axial middle of the pipe at the peak flow rate. The Lam- 
Bremhorst results follow the universal profile fairly closely. (The universal profile 
is based on steady flow.) Except at the axis, the Kim results do a good job of 
following the unsteady experimental data, at least for this instant in time. 

Figure 4.27 shows the variation of the friction coefficient over the first half cycle 
for the Lam-Bremhorst and Kim models. The experimental results of Friedman 
and the steady flow results at the corresponding bulk flow rate are also shown. The 
Kim results follow the laminar-like behavior during the first 30° crank angle, but 
then become turbulent well before the experimental data and about 10° before the 
Lam-Bremhorst results. Moreover, the Kim results are significantly higher than all 
others through most of the remainder of the half cycle. The friction coefficient at 
90° (peak flow) is very close to the experimental result, which explains the good 
velocity profile agreement shown by the Kim model in figure 4.26. Recall that 
the experimental velocity profile at 90° crank angle (figure 4.26) was influenced 
by the arrival of the turbulent slug. No traveling slug is expected when using the 
expansion/contraction geometry, so the close agreement between the Kim model 
and the experimental velocity profile results in figure 4.26 is fortuitous. 

The presence of a traveling slug being carried into the pipe after flow reversal is 
inferred from experimental observations within the pipe section itself. Turbulence 
levels are seen to rise dramatically at a time consistent with the arrival of a slug 
traveling through the pipe at the bulk fluid velocity. No detailed experimental data 
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is available concerning the structure of turbulence generated within the nozzle dur- 
ing outflow. On the other hand, the numerical models provide detailed information 
about the distribution of all variables within the nozzle. This information can be 
studied to observe the ability of the model to predict the traveling slug effect. The 
Kim model was used to produce figures 4.28 through 4.33, which show surface plots 
of the turbulent viscosity within the right-end nozzle for several steps in time, be- 
ginning at flow reversal. In these figures, the pipe ends at an axial length of 64. The 
nozzle is represented by the area to the right. In figure 4.28, the bulk flow has come 
to rest. The turbulence within the nozzle results from separation which occurred 
during the previous quarter-cycle, for which flow was to the right. As the cycle 
continues, this turbulence, characterized in these figures by the turbulent viscosity, 
is expected to be convected back into the pipe (to the left). Figure 4.29 shows the 
turbulent viscosity in the nozzle at 2.5° after flow reversal. The bulk flow is moving 
into the pipe at Re = 518. The levels of turbulent viscosity have decayed slightly, 
but there is no evidence that the turbulent structure has been convected from the 
nozzle into the pipe. Figures 4.30 through 4.33 show the turbulent viscosity at ap- 
proximately 5.7°, 9.2°, 13.0° and 17.1° crank angle after flow reversal, for which the 
instantaneous Reynolds numbers are 1187, 1926, 2709 and 3519, respectively. The 
turbulence in these figures is seen to decay in place, instead of being swept into the 
pipe. Similar flow visualizations based on the Lam-Bremhorst model (not shown) 
produced similar results. Both the Lam-Bremhorst and Kim models fail to ingest 
the turbulent slug after flow reversal. 

A further test was performed to determine whether the Kim model can predict 
the passage of a turbulent slug. The flow behind a turbulence generating grid was 
considered. If the turbulence level at the plane of the grid is increased and then 
decreased in a stepwise fashion, the turbulence within the fluid passing downstream 
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of the grid should exhibit a slug flow behavior. This type of flow is spatially one- 
dimensional. Near-wall effects need not be modeled and there is no turbulence 
generation, so the results axe a test of the convection and rate of dissipation of the 
turbulence. This test was also performed using the k-e model for comparison. 

Figure 4.34 shows the trajectory of the normalized turbulent kinetic energy 
downstream of a turbulence generating grid as produced by the Kim model. The 
length and time scales axe also normalized using k?/e and k/e, respectively. The 
inlet turbulence intensity is 2% initially, and this steady flow result is represented 
by the r = 0.0 curve. The inlet turbulence is then increased to 10% for the next 
five time steps before returning to the initial level. A slug of highly turbulent fluid 
is seen passing to the right after the inlet turbulence intensity is reduced, but this 
turbulence quickly dissipates while being convected downstream. 

Figure 4.35 shows the same results for the k-e model. For this model the turbu- 
lent slug persists slightly longer than in the case of the Kim model. For instance, 
the Kim model predicts a peak value of the normalized turbulent kinetic energy of 
about .0024 when r = 2.48, the first time step after the inlet turbulence intensity 
is reduced. At the same instant, the standard k-e model predicts a peak turbulent 
kinetic energy of about .0027. It appears that two time scale model does not provide 
any possibility of improvement over the k-e model with respect to the convection of 
a slug of turbulent fluid. 

The RNG-based k-e turbulence model was tested next for application to turbu- 
lent and transitional oscillating flow. The details of the High Reynolds Number 
(HRN) version of this model were discussed in chapter 2. 
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4.2.4 Results of the RNG Model 


The RNG-based k-e model of Yakhot and Smith (1992) has been tested to determine 
its suitability for the prediction of oscillating flow. The model was first tested 
on steady pipe flow using the HRN formulation. This model was not applied to 
oscillating flow due to the limitations of the wall function boundary conditions 
which were discussed in chapter 2. Since the details of the LRN form of the RNG 
model are not available in the literature, an attempt was made to formulate a low- 
Reynolds number RNG model by following the treatment of the Lam-Bremhorst 
model. The HRN RNG model was modified by applying the Lam-Bremhorst low 
Reynolds number functions to the turbulent viscosity equation (/ M ) and the first 
two source/sink terms of the dissipation equation (/i and /j). These functions are 
given by equations (2.25) through (2.27). The turbulent viscosity is then given by 

*-/>/,' c„7 (4-12> 

The dissipation equation becomes 
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This new model was first apphed to steady pipe flow for testing. The model 
tended to produce a gradual decay of turbulence, ultimately converging to a laminar- 
like solution. Since no damping function was applied to the last term, it was sus- 
pected to be the cause of the poor performance. This term has the ability to change 
sign, depending on the value of 77 , so further investigation was in order. For the 
purposes of this discussion, the last term will be referred to as the RNG coefficient 
multiplying The RNG coefficient is then 
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To better understand the behavior of this term, its variation with tj was first 
plotted as shown in figure 4.36. The RNG coefficient is negative for r, less than 4.38 
( 7 / 0 ). This produces amplification of turbulence. For tj greater than 4.38, the RNG 
coefficient causes damping of turbulence. To determine the variation of R and tj 
across the pipe section, their values were processed from the results of a converged 
data set using the Lam-Bremhorst model. The data set is based on steady pipe 
flow at Re = 12000, and the results are taken at the end of a 100 diameter long 
pipe. Figure 4.37 shows the profile of tj. The wall boundary condition for k forces 
tj to zero at the wall. The modest turbulence activity near the axis produces a 
small value of rj in that area. The value of tj is highest in the turbulence-generation 
region. Figure 4.38 shows the profile of R across the pipe section. Clearly, R will 
tend to enhance turbulence over the bulk of the pipe cross section. However, near 
the wall the effect of R is to cause significant damping. The magnitude of R in 
this region is such that it dominates the dissipation equation, forcing the solution 
toward a laminar-like result. This behavior must be modified in order to produce a 
plausible LRN RNG model. 

Several attempts were made to apply a damping function such as f 2 to R, but 
these efforts did not produce a well-behaved model. Finally, it was found that if the 
value of 7 , were restricted, R would be bounded and the model would produce the 
desired level of turbulence. The treatment adopted here is to restrict tj to values less 
than or equal to 4.0. This prevents R from contributing any turbulence damping 
behavior. While admittedly ad hoc, the model does produce good results for this 
steady pipe flow application. These results will now be discussed. 

Figure 4.39 compares the axial velocity profiles of the LRN RNG and Lam- 
Bremhorst models for steady pipe flow at Re = 12000. The RNG model produces 
a slightly higher peak and a rounder shoulder. The Lam-Bremhorst model is flatter 
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due to a somewhat overstated turbulent viscosity near the axis. The turbulent 
viscosity is compared in figure 4.40 for the same operating conditions. The reduction 
in turbulent viscosity at the axis for the RNG model is caused by the coefficient R 
in the dissipation equation. Near the axis, this coefficient acts as a sink term in the 
dissipation equation. Since R goes to zero at the axis, dissipation is boosted. This is 
an improvement over the Lam-Bremhorst model. Since turbulence generation goes 
to zero at the axis, turbulent viscosity is expected to drop, and this is reflected in the 
data of Reichardt (1951). Figure 4.41 compares the kinetic energy profiles. These 
are essentially identical near the axis, with the Lam-Bremhorst model producing 
a slightly higher peak. Near the wall, however, the RNG model produces slightly 
higher values. Figure 4.42 shows the dissipation profiles. Once again the Lam- 
Bremhorst model produces a slightly higher peak, but the RNG model is somewhat 
higher near the axis, which produces the desired drop in turbulent viscosity. Figure 
4.43 compares the velocity profiles in wall coordinates. The agreement between the 
models is good, with the RNG model producing a slightly higher “wake near the 

axis. 

Figures 4.44 through 4.47 compare the axial velocity, turbulent viscosity, kinetic 
energy and dissipation profiles of the LRN RNG and Lam-Bremhorst models for 
steady pipe flow at Re = 100000. The axial velocity profile of the Lam-Bremhorst 
model is even flatter relative to the LRN RNG model as compared to the Re = 
12000 case as shown in figure 4.44. This is due to the increasing spread between the 
turbulent viscosities of the two models. For Re = 12000 the Lam-Bremhorst model 
produces a turbulent viscosity at the axis which is about 35% higher than the LRN 
RNG model. At Re = 100000, the spread is about 39%, as shown in figure 4.45. 
The kinetic energy profiles of the two models are quite similar at Re = 100000, as 
shown in figure 4.46. The LRN RNG model produces slightly smaller values near 
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the axis and the peak, and somewhat larger values very close to the wall, but the 
agreement is better than in the Re = 12000 case. For Re = 100000 the peak in the 
dissipation is very high and close to the wall, as shown in figure 4.47, from which 
it is difficult to distinguish between the LRN RNG and Lam-Bremhorst models. 
Figure 4.48 shows the comparison in the vicinity of the wall. Once again, the Lam- 
Bremhorst model produces a higher peak. The wall value of the dissipation is also 
higher in the Lam-Bremhorst model. Figure 4.49 compares the velocity profiles in 
wall coordinates. For Re = 100000 both models follow the universal profile closely. 
Again, the RNG model produces a slightly higher “wake” near the axis. 

The LRN RNG model presented here does a good job reproducing the details of 
steady pipe flow. It remains to be seen whether the model offers any improvement 
in the prediction of oscillating pipe flow. As a test, the model was applied to a 100 
diameter long straight pipe subjected to oscillating flow with Re max = 12000 and 
Va = 80, the SPRE operating conditions. The model was reasonably well behaved, 
but produced results which are strikingly similar to those of the Lam-Bremhorst 
model. This may have been anticipated, given that the LRN form of the RNG model 
borrowed heavily from the Lam-Bremhorst model. Figure 4.50 shows the friction 
coefficient at the axial middle of the pipe over the first half cycle. The experimental 
and steady flow results are shown for reference. The LRN RNG results are virtually 
identical to those of the Lam-Bremhorst model (see figure 4.27). 

The LRN RNG k-e model tested above is not representative of the LRN version 
proposed by Yakhot and Orszag (1986). That model employed differential relations 
for the turbulent viscosity in the near-wall region. Details of the application of the 
model have not been published. Testing of the LRN RNG model of Yakhot and 
Orszag will have to wait until such time as this model becomes available. 
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Figure 4.4: Velocity profiles at the pipe midlength at select times during the flow 
cycle for V a = 60. 
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Figure 4.5: Velocity profiles at the pipe midlength at select times during the flow 
cycle for Va = 80. 
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Figure 4.6: Velocity profiles at the pipe midlength at select times during the flow 
cycle for V a — 100. 
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Figure 4.13: Friction coefficient over one half flow cycle at the pipe midlength for 
Va = 80. 
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Figure 4.14: Friction coefficient over one half flow cycle at the pipe midlength for 
Va = 100. 


91 



Value 



Crank Angle 

A = cos[crank ang] B = cosjcrank ang + 180] 


Figure 4.15: Variation of an acceleration source term for the dissipation equation, 
proposed by Koehler (1990). 
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Crank Angle = 30.0 Crank Angle = 210.0 



Crank Angle = 121.3 Crank Angle = 301.3 



Crank Angle = 180.0 Crank Angle = 360.0 

Figure 4.16: Streamlines within the left end region at select times during the flow 
cycle for V a = 40. 
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Crank Angle = 30.0 Crank Angle — 210.0 
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Crank Angle = 180.0 Crank Angle = 360.0 


Figure 4.17: Streamlines within the left end region at select times during the flow 
cycle for Va = 60. 
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Crank Angle = 121.3 Crank Angle = 301.3 



Crank Angle = 180.0 Crank Angle = 360.0 


Figure 4.18: Streamlines within the left end region at select times during the flow 
cycle for V a — 80. 
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Crank Angle = 180.0 Crank Angle = 360.0 


Figure 4.19: Streamlines within the left end region at select times during the flow 
cycle for V a = 100. 
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Figure 4.20: Grid for nozzle-ended model. 
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Figure 4.21: Friction coefficient during the first half of the flow cycle for the nozzle- 
ended model. 
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Figure 4.24: Normalized turbulent kinetic energy profiles 2 k/U^ max obtained from 
the Lam-Bremhorst and Kim turbulence models. 
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Figure 4.26: Axial velocity profiles in wall coordinates for the Lam-Bremhorst and 
Kim turbulence models. Experimental data and universal profile shown for refer- 
ence. 
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Figure 4.27: Friction coefficient over one half flow cycle at the pipe midlength for 
the Lam-Bremhorst and Kim turbulence models. Experimental data and steady 
flow results shown for reference. 
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Figure 4.28: Turbulent viscosity within smooth nozzle end at flow reversal. 
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Figure 4.29: Turbulent viscosity within smooth nozzle end 2.5° after flow reversal. 
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Figure 4.30: Turbulent viscosity within smooth nozzle end 5.7° after flow reversal. 
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Figure 4.31: Turbulent viscosity within smooth nozzle end 9.2° after flow reversal. 
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Figure 4.32: Turbulent viscosity within smooth nozzle end 13.0° after flow reversal. 
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Figure 4.33: Turbulent viscosity within smooth nozzle end 17.1° after flow reversal. 
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Figure 4.35: Normalized turbulent kinetic energy behind turbulence generating grid 
for Lam-Bremhorst model. 
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Figure 4.37: Profile of rj for Re = 12000, based on the Lam-Bremhorst data set. 
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Figure 4.39: Comparison of the axial velocity profiles of the RNG and Lam- 

Bremhorst models for Re = 12000. 






Figure 4.41: Compaxison of the kinetic energy profiles of the RNG and Lam 
Bremhorst models for Re = 12000. 
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Figure 4.44: Comparison of the axial velocity profiles of the RNG and Lam 
Bremhorst models for Re = 100000. 
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Figure 4.46: Comparison of the kinetic energy profiles of the RNG and Lam 
Bremhorst models for Re = 100000. 
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Figure 4.50: Friction coefficient over one half flow cycle at the pipe midlength for 
RNG model with Re max = 12000, Va = 80. 


127 


Chapter 5 

The Laminar Heat Transfer 
Solution 


In this chapter the heat transfer results for laminar pipe flow will be presented. 
Fully developed results will be considered first 3 s a means of verifying the computer 
model. The unsteady case will then be studied to determine the effect of flow 
oscillation on the heat transfer. 

5.1 Fully Developed Laminar Heat Transfer 

Fully developed laminar flow was considered in Chapter 3. In the case of constant 
property flow, the corresponding heat transfer problem can be solved conveniently 
once the flow solution is available. Two types of boundary conditions are commonly 
encountered in analytical heat transfer solutions: the fixed wall temperature and 
the prescribed wall flux. These limiting cases can be approximated in the laboratory 
by exposing the outer pipe wall to condensing steam or by wrapping the pipe with 
resistant heating wire. Most engineering systems may be approximated using one or 
the other method of heating, and in practice most situations lie somewhere between 
these limiting cases. 

The fluid temperature will necessarily change in the axial direction when heat 
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transfer occurs between the pipe wall and the fluid. However, a thermally fully 
developed situation is said to occur when the following condition is satisfied 


d_ 

dx 


[ T w (x) ~r(r,x) j _ (5.1) 

T w (x) - T b (x) 

This condition may occur when the boundary condition is either a fixed wall tem- 
perature or a prescribed wall flux. Boundary conditions for thermally developed 
flows are described in detail in Sparrow and Patankax (1977). 

The thermally fully developed laminar heat transfer problem can be solved an- 
alytically to determine the heat flow from the pipe inner wall to the fluid. When 
the local Nusselt number is defined as 


j.t hD 
Nu = —r~ 

k 


Incropera and DeWitt (1985) give the numerical values of 


Nu = 3.66 


(5.2) 


(5.3) 


for the fixed wall temperature boundary condition, and 


Nu = 4.36 


(5.4) 


for the prescribed flux boundary condition. 

The computer model was tested for heat transfer in laminar pipe flow to en- 
sure that these values were produced. Figures 5.1 and 5.2 show the evolution of 
the Nusselt numbers for the fixed wall temperature and prescribed flux boundary 
conditions, respectively, for developing flow. Each Nusselt number has been nor- 
malized by the corresponding thermally fully developed value given above. The 
flow is uniform across the inlet plane, so both the flow and heat transfer solutions 
require some distance before becoming fully developed. Figures 5.1 and 5.2 show 
that the predicted Nusselt numbers do approach the analytical values after a cer- 
tain development length. The hydrodynamic development length given by equation 
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(3.1) when Re = 2000 is 100 pipe diameters. The heat transfer solution should be- 
come thermally fully developed at some point after the flow solution. The computer 
predicted Nusselt numbers are very nearly fully developed after 100 pipe diameters, 
and axe essentially equal to the fully developed values after 200 diameters. This 
serves to verify the heat transfer portion of the computer program. 

5.2 Heat Transfer in Laminar Oscillating Pipe 
Flow 

The fluid flow solution for laminar oscillating pipe flow was considered in chapter 
3. In this section we consider the corresponding heat transfer solution. The heat 
transfer solution may be obtained by solving the energy equation (2.4) using the 
velocity field as an input. 

The computer program was configured to model the heat transfer by setting 
the pipe wall and expansion /contraction shoulders to a fixed hot temperature. The 
outer diameter of the expansion/contraction regions is adiabatic. For the laminar 
runs the pipe was 50 diameters long and the expansion/contraction regions were 
each 5 diameters long. In the SPRE Stirling engine the heater pipe outer walls 
are bathed in a liquid metal reactor coolant. The metal pipe walls will have a 
thermal conductivity much higher than the working fluid. For these reasons the 
fixed temperature wall boundary condition was selected as a good approximation. 
The inlet fluid was assigned a cool temperature. Heat was thus transferred from 
the pipe and shoulder walls to the fluid, much as would be expected in a shell and 
tube heat exchanger. 

The dimensionless temperature profiles at the axial middle of the pipe at select 
crank angles during the first half of the flow cycle are shown in figures 5.3 through 
5.5 for V a = 1, 30 and 100. Figure 5.3 shows that the unsteady effect is slight for 
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this low oscillation rate. The temperatures at the middle of the pipe are highest 
when the flow rate is low and the residence time of the fluid within the pipe is 
high. As the flow increases the temperatures drop. At 90° crank angle the fluid 
near the axis has had little opportunity to heat up before reaching the middle of 
the pipe. During the second quarter of the flow cycle the fluid decelerates and the 
temperatures rise once again. The effect of oscillation is apparent in the difference 
in temperature levels during accelerating and decelerating phases of the flow cycle. 
Early in the first quarter of the flow cycle the fluid has been heated significantly 
during the previous flow reversal. In the next quarter cycle the fluid is cooler since 
the pipe has been flushed with cool incoming fluid flowing at relatively high speed 
through the pipe. If the oscillation rate were very low (Va < 1) the profiles for 
crank angles of 30.0° and 150° would be identical, as would the profiles for 58.7° 
and 121.3°. 

The unsteady effects become more complicated when V"a is increased to 30, as 
shown in figure 5.4. For this higher rate of oscillation the fluid near the middle of 
the pipe never reaches the ends during a complete flow cycle. In other words, cool 
fluid is never swept entirely through the pipe. The temperatures near the middle 
of the pipe are thus noticeably higher than in the Va = 1 case. The variation of 
temperature throughout the cycle is so complex that a lucid commentary is difficult 
to provide. This is caused by the complicated velocity profiles seen previously in 
figure 3.9. The variation of the velocity near the wall is ahead of, and near the axis 
is behind, that of the bulk flow. 

The temperature profiles become very flat when Va is raised to 100, as shown 
in figure 5.5. The fluid particles at the axial middle of the pipe travel only a small 
axial distance during the course of a flow cycle. 

The temperature contours within the pipe at select times during the first half 
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cycle are shown in figures 5.6 through 5.8 for Va — 1, 30 and 100. The radial 
dimension of each subplot has been stretched by a factor of 8 for clarity. The top of 
each subplot represents the pipe wall, where heat is introduced, while the bottom 
represents the pipe axis. The flow is from left to right in all subplots except the last, 
which represents bulk flow reversal. The temperature contours during the second 
half of the cycle are mirror images of those shown for the first. Figure 5.6 reveals 
temperature contours which are quasi-steady in nature. The cool incoming fluid is 
heated progressively as it passes through the pipe. The speed of the flow can be 
gauged by the length to which the cool incoming fluid penetrates into the pipe. At 
peak flow, for instance (90° crank angle), the lowest temperature contour meets the 
axis nearly three-quarters of the way through the pipe. For all other times the fluid 
is heated much sooner. 

The temperature distribution is more complicated when the oscillation rate in- 
creases. Figure 5.7 shows the contours for Va = 30. Two interesting features deserve 
attention in the first subplot of figure 5.7. First is the complicated distribution at 
the left end of the pipe for 29° crank angle. This is caused by the recirculation which 
was set up in the left end region during the previous cycle. At that time, a bubble 
of heated fluid was trapped in the end region. When the present cycle started, the 
slow moving incoming fluid was forced to climb around the bubble. Compare with 
figure 3.13. As a result, the coolest fluid enters the pipe at a radius between the 
wall and the axis. The second feature is the presence of cooler fluid at the right end 
of the pipe. This fluid was swept into the pipe during the end of the last cycle, and 
has had less time to be heated relative to the fluid near (but not at) the left end. 
This effect becomes more pronounced as the oscillation rate increases. 

The temperature contours for Va = 100 are shown in figure 5.8. The situation 
is qualitatively similar to the Va = 30 case, but the fluid is much warmer. Also, the 
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warmest fluid along the pipe length is located near the axial middle. In fact, this 
very warm fluid never reaches the end of the pipe, a feature which is characterized 
by the amplitude ratio of the bulk fluid motion A r . At this high rate of oscillation 
the amplitude ratio is only 0.2, showing that a fluid particle moving at the bulk 
fluid velocity would be displaced in the axial direction a distance of only two-tenths 

of the pipe length during the course of a flow cycle. 

The cycle-averaged variation of the dimensionless bulk temperature is shown in 
figure 5.9 for Va = 1, 30 and 100. This type of figure is suggestive of the axial 
diffusion of heat which is driven by the axial temperature gradient. Figure 5.9 
shows that the axial bulk temperature gradient is modest for Va = 1, and thus the 
axial heat transfer is largely by convection. The temperatures at the pipe ends are 
higher than T tn since the expansion and contraction regions hold a volume of heated 
fluid. The relative effectiveness of diffusion increases with the rate of oscillation. 
The axial temperature gradient is larger when Va is increased to 30, as required m 
order to move heat along the pipe axis in the face of diminished convection. This 
effect is even more pronounced when Vo is increased to 100 and A r is only 0.2. The 
contribution of convective heat transfer is small under these conditions, so the axial 
transport of heat becomes dominated by gradient diffusion. 

The Nusselt number variation along the pipe wall at select crank angles during 
the first half of the flow cycle is shown in figures 5.10 through 5.12 for Va = 1, 30 
and 100. The flow is from left to right. Each curve is normalized by the Nusselt 
number for fully developed laminar heat transfer with a constant wall temperature 
boundary condition (Nu fd = 3.66). For low Va , the flow and heat transfer problem 
is quasi-steady. The variation of the Nusselt numbers shown in figure 5.10 is similar 
to that of figure 5.1. The developing temperature field near the pipe entrance 
creates a Nusselt number well above the fully developed value. The Nusselt number 
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then drops and levels off as the end effect diminishes. The hydrodynamic (and 
heat transfer) development length increases with Re, so the level of Nu depends 
slightly on the crank angle (bulk flow rate). The only unusual feature of figure 5.10 
is the fact that the curves of Nu are lower during the first quarter of the cycle. For 
instance, the curve for 30° crank angle is lower than for 150°. This hysteresis is 
caused by the oscillating flow effect. The wall heat flux is low after flow reversal 
due to the extensive heating of the near-wall fluid which occurred while the flow 
rate was low, during the previous flow reversal. During the second quarter of the 
cycle the near-wall temperature is lower as a result of the significant convection of 
cool incoming fluid. This effect will also be seen in figure 5.16 below. 

Oscillation has a more complicated effect on the relative level of Nu during 
the cycle as Va increases. For crank angles early in the cycle, the near-wall flow 
responds quickly to the imposed pressure gradient, creating a laxge velocity gradient, 
as shown previously in figure 3.9 for a moderate Vo of 30. This in turn increases the 
temperature gradient and Nu. This effect is quite pronounced for xj D greater than 
about 15 when V a — 30, as shown in figure 5.11. Later in the cycle the gradients 
at the wall are lower. For Va = 30, flow reversal at the wall occurs at a crank angle 
near 149°, and the heat transfer is lowest at that time. 

When Va becomes large, thermal convection diminishes in comparison to dif- 
fusion. The bulk fluid motion is restricted, and A r is low. For laxge Va , a fluid 
particle near the middle of the pipe will never be swept out of the pipe, but will 
instead shuttle back and forth over a short distance. As a result, this fluid will reach 
a temperature closer to that of the wall. The heat flux and Nu are thus lowest near 
the middle of the pipe. This effect is shown in figure 5.12 for Va = 100. 

The variation of the axial-averaged Nu (Nu ) over the entire cycle is shown in 
figures 5.13 through 5.15 for Va = 1, 30 and 100. There are two cycles of heat 
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transfer in each flow cycle since the heat transfer does not depend on the axial 
direction of flow. The Nusselt number is normalized against Nuf d . Figure 5.13 
shows the time variation of Nil when Va is low. When the flow rate is large, the 
entire length of the pipe falls within the developing region and Nu is higher than 
Nu/d . The developing region takes up a smaller fraction of the pipe when the flow 
rate is low, and Nu approaches Nu Jd . The rapid variation near bulk flow reversal 
is an unsteady effect caused by circulation of fluid within the pipe. This effect 
diminishes shortly after flow reversal (20° and 200° crank angle) as the bulk flow 
increases and once again dominates the heat transfer. 

The circulation effect becomes more prominent as Va increases. Figure 5.14 
reflects the fact that near-wall flow reversal begins earlier and is significant over a 
larger portion of the cycle when Va = 30. This effect continues to increase with Va 
as shown in figure 5.15, for which Va = 100. 

The Nusselt number used thus far is defined as 



The heat transfer coefficient was defined as 


(5.5) 


h = 


g 

T w — Tb 


(5.6) 


where q is the local heat flux from the wall to the fluid, T w is the wall temperature, 
and Tj, is the local bulk fluid temperature. The term local refers to a quantity whose 
value depends on axial position. 

This definition of Nu is common for steady flow heat exchanger design, for which 
Tb is used in place of the free-stream temperature (Too) used for external flows. 

The bulk (velocity- area weighted) temperature for incompressible, constant prop- 
erty flow has been defined here as 

_ / |u(r)|Trdr ( 5 

b /|u(r)|rdr 
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where r dr is a differential flow area and u(r) is the axial velocity component. The 
absolute value of velocity is used since the scalar Tj, is sensitive only to the magnitude 
of the flow. This definition also applies at flow reversal. Without the use of the 
absolute value the Tb definition would produce Nu = 0 when the bulk velocity 
becomes zero. 

Figures 5.16 through 5.18 show the phase relationship between flux, T w - Tb , 
and Nu over one flow cycle at the axial middle of the pipe for V a = 1, 30 and 
100. This type of comparison serves to further illuminate the complicated features 
of this unsteady flow and heat transfer problem. The sudden jump in Nu near flow 
reversal has been discussed above. The underlying variation of flux and T w - Tb are 
seen clearly in figure 5.16. As flow reversed approaches, the temperature difference 
drops faster than the flux, boosting Nu . At this low oscillation rate the flow is 
quasi-steady, but the flux and temperature difference do not have the same trends 
throughout the cycle. Except for a short period following flow reversal, the flux 
tends to rise later than the temperature difference during accelerating portions of 
the cycle and fall sooner during decelerating portions. This behavior was discussed 
above in conjunction with figure 5.10 and tends to keep Nu slightly lower during 
acceleration. Note that the variation of the flux and temperature difference at low 
oscillation rate bear some resemblance to the absolute value of a sinusoid. The flux 
and temperature difference both peak at nearly the same instant as the flow, and 
reach their mi nim a together only shortly after the flow. 

The variation of the flux and temperature difference become more complicated 
as Va increases. Figure 5.17 shows that the flux and temperature difference lag 
behind the flow, peaking at about 156° and 165°, respectively, for Va = 30. The 
heat flux variation is vaguely sinusoidal, though it rises faster during deceleration 
than it falls during acceleration. The variation of the temperature difference is 
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particularly complex due to the combination of phase shift and flow reversal effects. 
The situation becomes more complicated yet when Va = 100, as shown in figure 
5.18. Here the heat flux lags the flow by 90° and the variation appears nearly 
sinusoidal. 

Figures 5.16 through 5.18 consider the heat transfer behavior only at the axial 
middle of the pipe. At this point the end effects are small when the flow is low 
and symmetric in any case. The heat transfer behavior is even more complicated 
when viewed at positions away from the middle of the pipe. Figure 5.19 shows the 
variation at a position 10.8 diameters from the left end of the pipe for Vq 1, 
the same operating condition used in figure 5.16. At this position the heat transfer 
is significant during the first half of the cycle, when the left end of the pipe is 
the entrance. During the second half the left end is closer to the outlet, so the 
heat transfer is closer to the fully developed case. The differences in heat transfer 
behavior between the first and second halves of the flow cycle only become larger 
as either end of the pipe is approached. 

One of the objectives of the present research is to reduce the unsteady 2-D heat 
transfer behavior in oscillating pipe flow into a more manageable form that can be 
incorporated into spatially 1-D Stirling engine performance codes. The temporal 
and spatial variation of Nu seen in figures 5.10 through 5.19 is far too complicated 
to reduce to a form such as 

Nu = Nu(x,t) ( 5 - 8 ) 

It would be desirable to devise some means of simplifying the Nu expression so 
that its behavior can be expressed by a simple correlation. The unsteady and flow 
reversal effects discussed above produce a complicated variation of Tt , and this is 
the most ill-behaved part of the Nu definition. The use of Tt may not even be 
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desirable from the perspective of the engine performance modeler. The detailed 
variation of 7t through space and time will not be known a priori, but would be 
required in order to compute the total heat transfer between the fluid and pipe wall. 
The inlet and outlet fluid temperatures will be available from the performance code 
to be used (iteratively) as an input to the heat transfer calculation. It is proposed 
instead that the temperature difference T w - 7i„ be used in the definition of the 
heat tr ans fer coefficient given by equation (5.6). This will produce a variation of 
Nu which is identical to that of the flux, since T w - T, n is a constant. The resulting 
redefined Nusselt number will be referred to as Nu*. 

Removing the Tj, dependence eliminates some of the complications due to the 
oscillating flow effect. The end effects remain to be dealt with. We have seen that 
the end effects axe significant, but they can be accounted for implicitly by the use 
of axial averaging. The final correlation for Nu would then be restricted to a given 
L/D ratio but would still be a function of time. This is acceptable for the purpose 
of building a Stirling engine performance code once the heat exchanger tube has 
been selected. (Cost, weight, configuration and space constraints may be just as 
important as the heat transfer performance.) When the redefined Nusselt number 
is averaged over the length of the pipe (NtT) we find the behavior shown in figures 
5.20 through 5.22. 

Figure 5.20 shows the predicted variation of Nu over one flow cycle for Va — 
1. The prediction is shown in square symbols. A curve fit is also plotted using the 
equation 

~NxT = 4.0 * |sin(cranfc angle — 5°)| + 0.5 (5-9) 

Both the numerical prediction and the curve fit have a cycle- averaged Nu of 3.04. 
The curve fit was established by adjusting the amplitude, phase shift (relative to 


138 



the bulk flow) and vertical offset to produce the same cycle- averaged value and the 
best qualitative agreement with the computed curve. The form of the curve fit (5.9) 
is acceptable whenever Va is less than 5. 

Figure 5.21 shows the variation of TVtTover one flow cycle for Va = 30. For this 
case the following sinusoidal curve fit was selected 

Nu~ = 0.32 * sin (2 * ( crank angle — 90°)) + 2.24 (5.10) 

The cycle- averaged value is 2.24. 

Figure 5.22 shows the variation of Nu ~ over one flow cycle for Va = 100. The 
curve fit is given by 

~Nv* = 0.13 * sin (2 * ( crank angle — 69°)) + 1.43 (5.11) 

The cycle- averaged value is 1.43. 

Figures 5.20 through 5.22 show that, unlike the conventionally defined Nusselt 
number, Nu~ decreases with increasing Va . This clearly illustrates the reduced 
effectiveness of axial convective transport as the oscillation rate increases. At very 
high rates of oscillation the axial transport is primarily due to diffusion which 
has been enhanced through interaction with the oscillating flow field. Zhang and 
Kurzweg (1991) refer to this condition as enhanced axial heat transfer. The axial 
heat flux under these conditions may be substantially larger than in the case of pure 
conduction, but will nevertheless be smaller than if steady pipe flow were allowed 
to occur. 

The variation of the amplitude, mean and phase shift (relative to the bulk flow) 
of N u* with Va is shown in figure 5.23. This figure may be used as a means of 
interpolating the laminar Nu ’ correlation (equations (5.10) or (5.11)) for oscillation 
rates in the range 5 < Va < 100. 
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Figure 5.1: Development of the Nusselt number along the pipe length for a fixed 
wall temperature boundaxy condition. 
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Figure 5.2: Development of the Nusselt number along the pipe length for a pre- 
scribed flux boundary condition. 
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Figure 5.6: Temperature contours within pipe during first half cycle for Va = 1. 
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Figure 5.7: Temperature contours within pipe during first half cycle for Va = 30. 
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Figure 5.8: Temperature contours within pipe during first half cycle for Va = 100. 
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Figure 5.9: Cycle-averaged dimensionless bulk temperature for Va = 1, 30 and 100. 
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Figure 5.13: Axial-averaged Nusselt number variation during the cycle for Va — 1. 
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Figure 5.14: Axial-averaged Nusselt number variation during the cycle for V a = 30. 
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Figure 5.15: Axial-averaged Nusselt number variation during the cycle for V a = 

100 . 
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Figure 5.17: Phase compaxison between heat transfer quantities at the axial middle 
of the pipe for Va = 30. 
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Figure 5.23: Amplitude, mean and phase shift for the laminar Nusselt number 
correlation as a function of V a . 
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Chapter 6 

The Turbulent Heat Transfer 
Solution 


In this chapter the heat transfer results for turbulent flow will be presented. The 
steady results will be examined in order to verify the computer program. The 
uns teady results will then be studied to determine the effect of flow oscillation on 
the heat transfer in the presence of transition, turbulence and relaminarization. 

6.1 Fully Developed Turbulent Heat Transfer 

The axial variation of the Nusselt number in developing turbulent flow is shown in 
figure 6.1 for Reynolds numbers of 10000 and 50000. The calculations are based on 
a pipe with L/D = 100 and a uniform inlet velocity profile. The inlet turbulence 
intensity is 5%. Grid sizes of 32 by 42 and 32 by 52 were used for Re = 10000 and 
50000, respectively. The results are normalized against the Gnielinski correlation 
reported by Kakac, et al. (1987). This correlation is suitable for the ranges 2300 < 
Re < 5 * 10 6 and 0.5 < Pr < 2000, with errors on the order of 5% for modest 
Pr and 10% for larger Pr. The behavior of the two curves differs near the inlet of 
the pipe. For Re = 10000 Nu tends to drop below Nu Jd near the inlet. This is 
due to the specification of the turbulence kinetic energy and dissipation rate at the 
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inflow boundary, combined with the uniform inlet velocity profile and the modest 
Re. As the flow enters the pipe, the velocity gradient over the most of the pipe 
cross section is zero, resulting in zero turbulence generation. Turbulence is once 
again generated as the velocity profile develops. After a few diameters the level of 
Nu is within 0.7% of Nufd . This effect can also be observed when Re = 50000, but 
is barely visible in figure 6.1. The computed value of Nu for Re = 50000 is 4.0% 
high. It is not obvious whether the errors are due to the fully developed correlation 
or the computer model, but the differences are slight. The numerical results can 
be brought into agreement with the correlation by increasing the turbulent Prandtl 
number Pr t , which was assigned the uniform value of 0.9 for this calculation. This 
level of agreement can be considered as validation of the model for steady flow heat 
transfer. 


6.2 Heat Transfer in Turbulent Oscillating Pipe 
Flow 

The treatment used to study the turbulent oscillating heat transfer problem is simi- 
lar to that used in Chapter 5 for the laminar case, except that the Reynolds averaged 
energy equation (2.44) is used. A pipe is fitted with expansion and contraction end 
regions whose outer diameter is adiabatic. The pipe walls and the shoulders of the 
expansion and contraction end regions are set to a uniform hot temperature, while 
the inflow is set to a uniform cold temperature. The expansion and contraction 
regions are each 10 pipe diameters long. The Lj D ratio for the model is 60. See 
figure 2.1. The flow and oscillation rates are given by Re max = 12000 and Va = 80. 
These operating conditions match those of the heater section of the SPRE engine. 

The dimensionless temperature profiles at the axial middle of the pipe (x/D = 
30) at select crank angles during the first half of the flow cycle are shown in figures 
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6.2 through 6.5 for Va = 40, 60, 80 and 100. These figures show that the level 
of temperature at the axial middle of the pipe increases with V a . This is caused 
by the decreased convection which occurs as oscillation frequency is increased. The 
same effect was observed in figures 5.3 through 5.5 for laminar flow. The sequence of 
events that produces the temperature profiles shown in figure 6.2 can be described 
as follows. The fluid at the axial middle of the pipe at 30° crank angle has been in 
the pipe since prior to the last flow reversal and so has a fairly high temperature. 
The same is true at 58.7°, so the temperature rises further. By 90° crank angle 
fluid from outside the pipe has reached the axial middle of the pipe, causing the 
temperature to drop. The temperatures continue to drop until very nearly 180° 
crank angle, when convection is small and the temperatures within the pipe rise 
almost exponentially. For the given values of Re max , Va and L/D, the amplitude 
ratio A r for figure 6.2 is 2.5, showing that a significant amount of fluid is convected 
through the pipe during each half cycle. The amplitude ratio decreases in proportion 
as Va increases. In figure 6.3 A,, is reduced to 1.67 when Va is 60. As a result, cool 
fluid is only just beginning to arrive at the axial middle of the pipe by 90° crank 
angle. Similar conclusions can be drawn from figures 6.4 and 6.5 for Va = 80 and 
100 and A r = 1.25 and 1.0, respectively. 

The temperature contours within the pipe at select times during the first half 
cycle are shown in figures 6.6 through 6.9 for Va = 40, 60, 80 and 100. The radial 
dimension of each subplot has been stretched by a factor of 10 for clarity. The top 
of each subplot represents the pipe wall, where heat is introduced, while the bottom 
represents the pipe axis. The flow is from left to right in all subplots except the last, 
which represents bulk flow reversal. The temperature contours during the second 
half of the cycle are mirror images of those shown for the first. In figure 6.6 the 
temperature contours show a behavior similar to that of figure 5.7 for the laminar 
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case. Starting with the subplot for 30° crank angle, the progression of the fluid to 
the right can be seen. The warmest fluid is near the left end, having spent part 
of the last cycle being heated inside the pipe. As the cycle continues, this fluid is 
swept to the right, followed by cooler fluid. The very warmest slug of fluid which 
is seen in the pipe at 58.7° crank angle is swept out of the pipe prior to 90° crank 
angle. 

The temperature contours behave in a similar manner for V a = 60, 80 and 100, 
but the warm slug of fluid persists in the pipe for an increasing portion of the cycle. 
For the Va = 100 case, A r is 1.0. The fluid near the axis moves faster than the bulk 
fluid, but figure 6.9 shows that the warm slug only just has time to reach the end 
of the pipe before flow reversal in the Va = 100 case. The results for V a = 60 and 
80 are intermediate between those for Va = 40 and 100. 

The cycle-averaged dimensionless bulk temperature variation along the pipe is 
shown in figure 6.10 for Va = 40, 60, 80 and 100. The trend with increasing Va 
is the same as in the laminar case, shown in figure 5.9. As Va increases and the 
convective effect decreases, the fluid temperature must rise in order to move heat 
in the axial direction by means of diffusion. Figure 6.10 gives an indication of the 
net effectiveness of diffusion over a complete oscillating flow cycle. 

The Nusselt number along the pipe wall at select crank tingles during the first 
half cycle is shown in figures 6.11 through 6.14 for Va = 40, 60, 80 and 100. The 
flow is from left to right. Each curve is normalized by the Nusselt number for 
fully-developed turbulent heat transfer (Nujj) at the corresponding instantaneous 
Reynolds number (Re). The developing temperature field near the pipe entrance 
creates a Nusselt number which is well above the fully-developed value. The Nusselt 
number then drops and levels off as the end effect diminishes. The effect of oscil- 
lation appears in the change in the level of Nu during the cycle. For crank angles 
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early in the cycle, the flow retains a laminar-like character, with little turbulent 
cross-stream transport. As a result, Nu drops below Nu/d . As the flow increases, 
transition occurs and Nu rises. During deceleration, turbulence developed at peak 
flow persists and maintains a high cross-stream transport, boosting Nu above Nujd- 
This explanation needs further discussion for the case when Va = 100. In this case 
the rate of oscillation is sufficiently high that the flow retains some turbulent char- 
acter at 30° crank angle. This turbulence is still in the process of decaying from 
the peak levels which occurred during the last cycle. On the other hand, transition 
to turbulence in the present cycle is delayed by the acceleration effect until just 
after 60° crank angle, giving Nu a laminar-like behavior for 58.7°. This sequence of 
events is also observed in the behavior of the friction coefficient, shown in figures 
4.11 through 4.14. 

Figures 6.15 through 6.18 show the axial-averaged Nu over one flow cycle for 
Va = 40, 60, 80 and 100. There axe two cycles of heat transfer in each flow cycle 
since the heat transfer does not depend on the axial direction of flow. The Nusselt 
number is normalized against Nu/d at the corresponding Re. Nu is particularly 
large near flow reversal due to circulation caused by the oscillating flow effect. Nu 
drops below Nuj d early in the cycle due to the acceleration damping effect. Nu 
rises above Nujd later as discussed in conjunction with figures 6.11 through 6.14 
above. Notice that the minimum of each curve falls lower and occurs later as Va 
increases, once again due to acceleration damping of turbulence. 

The Nusselt numbers presented in figures 6.11 through 6.18 are based on the 
definition used in equation (5.5) 
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The heat transfer coefficient is given by equation (5.6) 


k T w -T b 

The bulk (velocity-area weighted) temperature for incompressible, constant prop- 
erty flow is defined by equation (5.7) 

£Kr)j Trdr 
b f |u(r)| r dr 

As was seen in the laminar cases, the wall heat flux q , the temperature difference 
T w - T b and the Nusselt number Nu vary in a complicated fashion over the course 
of the oscillating flow cycle. Figures 6.19 through 6.22 show the phase relationship 
between the flux, T w - T b , and Nu over one flow cycle at the axial middle of the 
pipe for Va = 40, 60, 80 and 100. Three features in these figures will be discussed 
further: 

1. the sudden change in behavior of T w - T b and Nu neax flow reversal 

2. the phase shift between the three curves 

3. the sharp rise in flux after each flow reversal. 

The behavior of T w - T b and Nu near flow reversal is an unsteady effect. As the 
bulk flow rate approaches zero near flow reversal, the convection has a diminished 
effect on the heat transfer problem. At flow reversal, the heat transfer is nearly 
reduced to a transient conduction problem, for which an exponential behavior is 
expected. If the flow were to remain at rest after completing a number of cycles, 
the value of T w - T b would exponentially approach zero as the fluid absorbed heat 
from the wall. Soon after flow reversal, the convection again becomes important, 
imparting a behavior which is periodic due to the prescribed bulk flow boundary 
condition. The value of Nu neax flow reversal also reflects this unsteady effect. 
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The peak in the curve T w - T b lags the peak flux by roughly 20° to 40° as Va 
increases from 40 to 100. This is also an unsteady effect. As the heat flux at the 
wall changes, T b cannot respond instantaneously, due to thermal inertia. The peak 
value of Nu , based on the ratio of flux to T w - T b , leads the flux by roughly 15° to 
50°. Moreover, the behavior of Nu is not sinusoidal. This effect can be seen in a 
test case which considers the ratio of two sinusoidal functions separated by a small 
phase difference. Figure 6.23 shows three functions A , B, and C, given by 


A = 2 + sin(6) 

(6.1) 

B — 2 + sin{6 <f > ) 

(6.2) 

C = A/B 

(6.3) 


where <f> gives the phase shift between A and B. In this test case, functions A, B 
and C may be considered analogues to the heat flux, T w - T b , and Nu , respectively, 
of figures 6.19 through 6.22. As shown in figure 6.23, function C leads both A and 
B and is not sinusoidal. The function C reaches its peak roughly 80° after each 
valley, with the next valley following 100° later. This phase shift effect also causes 
a non-sinusoidal Nu when Nu is based on a flux and temperature difference which 
are out of phase. The phase shift between Nu and flux can be eliminated by using 
T w - Ti n instead of T w - T b in the definition of Nu . This causes Nu to have the 
same behavior as flux since T w - T tn is a constant. 

The sudden rise in flux after flow reversal is due to transition to turbulence. This 
rise occurs from roughly 35° to 50° crank angle after flow reversal for Va ranging 
from 40 to 100. Figure 6.24 shows the turbulent kinetic energy profiles for select 
crank angles at the axial middle of the pipe when Va = 80. The turbulent kinetic 
energy rises sharply near the wall at about 45°, demonstrating that the rise in flux 
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near that crank angle is due to transition. The same transition also occurs at 225° 
in the second half of the flow cycle for V a = 80. 

As was seen in the laminar case in chapter 5, the temporal and spatial variation 
of Nu seen in figures 6.11 through 6.22 is far too complicated to reduce to a function 
of axial position and time as 

Nu = Nu(x,t) 

Once again an alternate definition of Nu will be used along with axial averaging in 
order to produce a more tractable expression for the oscillating heat transfer. The 
pipe inlet temperature T,„ will be used in place of 7j, in the temperature difference, 
relieving the engine modeler of the need to specify the complicated variation of 
the bulk temperature. The inlet and outlet fluid temperatures will be available 
from the performance code to be used (iteratively) as an input to the heat transfer 
calculation. The end effects can be accounted for implicitly by the use of axial 
averaging. The resulting correlation for Nu will be restricted to a given L/D ratio 
but will still be a function of time. This is acceptable for the purpose of building a 
Stirling engine performance code once the heat exchanger tube has been selected. 
(Cost, weight, configuration and space constraints may be just as important as the 
heat transfer performance.) When the redefined Nusselt number is averaged over 
the length of the pipe ( Nu ) we find the behavior shown in figures 6.25 through 
6.28. 

Figure 6.25 shows the predicted variation of Nu ' over one flow cycle for V a = 
40. The prediction is shown in square symbols. A curve fit is also plotted using the 
equation 

Nu = 10.0 * sin (2 * ( crank angle — 68°)) + 14.93 (6-4) 

Both the numerical prediction and the curve fit have a cycle-averaged jVtt’of 14.93. 
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The curve fit was established by adjusting the amplitude, phase shift and vertical 
offset to produce the same cycle-averaged value and the best qualitative agreement 
with the computed curve. 

Figure 6.26 shows the variation of Nu over one flow cycle for Va — 60. For this 
case the following sinusoidal curve fit was selected 


~Nu = 7.5 * sin (2 * ( crank angle — 78°)) + 12.58 (6-5) 

The cycle-averaged value is 12.58. 

Figure 6.27 shows the variation of Nu over one flow cycle for Va = 80. The 
curve fit is given by 

Nu" = 5.5 * sin (2 * (crank angle — 88°)) + 9.72 (6.6) 

The cycle-averaged value is 9.72. 

Figure 6.28 shows the variation of Nu over one flow cycle for Va = 100. The 
curve fit is given by 

Nu * = 4.1 * sin (2 * ( crank angle — 85°)) + 8.17 (6-7) 


The cycle- averaged value is 8.17. 

Figures 6.25 through 6.28 show that, unlike the more conventional Nusselt num- 
ber defined by equations (5.5) through (5.7), TVV decreases with increasing Va . 
This clearly illustrates the reduced effectiveness of axial convective transport as the 
oscillation rate increases. At very high rates of oscillation the axial transport is 
primarily due to diffusion which has been enhanced through interaction with the 
oscillating flow field. The axial heat flux under these conditions may be substan- 
tially larger than in the case of pure conduction, but will nevertheless be smaller 
than if steady pipe flow were allowed to occur. 
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The variation of the amplitude, mean and phase shift (relative to the bulk flow) 
of jVit' with Va is shown in figure 6.29. This figure may be used as a means of inter- 
polating the turbulent 7hT correlation (equations (6.4) through (6.7)) for oscillation 
rates in the range 40 < Va < 100. 
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Figure 6.6: Temperature contours within pipe during first half cycle for Va = 40. 
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Figure 6.7: Temperature contours within pipe during first half cycle for Va = 60. 
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Figure 6.8: Temperature contours within pipe during first half cycle for Va = 80. 
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Figure 6.9: Temperature 
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Figure 6.12: Nusselt number versus axial position during the first half cycle for V a 
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Figure 6.13: Nusselt number versus axial position during the first half cycle for V a 
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Figure 6.14: Nusselt number versus axial position during the first half cycle for V a 

= 100 . 
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Figure 6.15: Axial- averaged Nusselt number variation during the cycle for Va = 40. 
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Figure 6.16: Axial-averaged Nusselt number variation during the cycle for Va = 60. 
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Figure 6.17: Axial-averaged Nusselt number variation during the cycle for Va = 80. 


189 



u.u-1 1 — 

0.0 45.0 


— l 1 i i i 

90.0 135.0 180.0 225.0 270.0 

Crank Angle 


315.0 


360.0 


Figure 6.18: Axial-averaged Nusselt number variation during the cycle for Va = 

100 . 
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Figure 6.19: Phase comparison between heat transfer quantities at the axial middle 
of the pipe for V a = 40. 


191 


Normalized Value 



0.0 45.0 90.0 135.0 180.0 225.0 270.0 315.0 360.0 

Crank Angle 


Figure 6.20: Phase comparison between heat transfer quantities at the axial middle 
of the pipe for V a = 60. 
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Figure 6.21: Phase comparison between heat transfer quantities at the axial middle 
of the pipe for V a — 80. 
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Figure 6.22: Phase comparison between heat transfer quantities at the axial middle 
of the pipe for V a — 100. 
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Chapter 7 

Concluding Remarks 


7.1 Summary and Conclusions 

1. A literature survey shows that sufficient experimental data is available for 
validation of numerical flow and transition predictions. 

2. Only limited data is available for validation of heat transfer results. The data 
is not suitable for comparison with the present work due to differences in the 
end conditions. 

3. The Low Reynolds Number k- e turbulence model of Lam-Bremhorst (1981) 
predicts tr ans ition and relaminanzation of the flow, but not at the same times 
observed during experiments. Previously proposed modifications to this tur- 
bulence model were examined, and found to be physically incompatible with 
oscillating flows (flows with reversal). The model does not produce the exper- 
imentally observed convection of large turbulent slugs. 

4. The two time scale turbulence model of Kim (1992) produces results which 
are similar to the Lam-Bremhorst model. No improvement was found in terms 
of transition, relaminarization or the convection of turbulent slugs. 

5. The High Reynolds Number RNG-based k-e model of Yakhot and Smith 
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(1992) provides modest improvements over Lam-Bremhorst for steady pipe 
flow, but cannot be applied to oscillating flow due to limitations of the wall 
functions boundary treatment. An extension of the Yakhot and Smith model 
to low Reynolds numbers was devised, but produced results nearly indistin- 
guishable from Lam-Bremhorst, from whom the low Reynolds number func- 
tions were borrowed. The rigorously derived Low Reynolds Number RNG 
model of Yakhot and Orszag (1986) is not yet fully described in the literature, 

and was not tested. 

6. The heat transfer predictions show that the wall heat flux and bulk tem- 
perature are out of phase with the bulk flow rate (and each other). The 
temperature difference exhibits complicated behavior near bulk flow reversal 
due to the momentary “conduction-like” nature of the heat transfer problem 
in the absence of (significant) convection. 

7. A Nusselt number based on the flux and bulk temperature is even less well 
behaved than the bulk temperature, with non-sinusoidal variation throughout 
the cycle and strong dependence on end effects. 

8. A newly defined Nusselt number based on the pipe inlet temperature instead 
of the bulk temperature is proposed. Also, axial averaging of this Nusselt 
number produces a more tractable expression for the oscillating heat transfer. 

7.2 Contributions of This Work 

This work is the first of which the author is aware which provides an investigation 
and results for the following: 
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1. Detailed numerical predictions for turbulent and transitional oscillating pipe 
flow with complicated end geometries. 

2. Application of the Kim two time scale turbulence model to highly unsteady 
flow. 

3. N um erical predictions for the variation of the Nusseit number along the pipe 
and throughout the oscillating flow cycle. 

4. A proposed modification to the Nusseit number, producing a form which is 
more readily applied to Stirling engine design. 

7.3 Suggestions for Further Research 

The following items deserve further attention during future research: 

1. When available, the Low Reynolds Number version of the RNG-based k-t 
model should be applied to the oscillating flow and heat transfer problem to 
determine whether it yields improvements over the Lam-Bremhorst model. 

2. The parametric range of this study should be expanded to ensure the results 
axe applicable to a broader range of oscillating flow and heat transfer problems. 

3. The geometry and boundary conditions should be adjusted to produce results 
which may be compared against the newly available experimental heat transfer 
data of Simon and Qiu (1993). 

4. The variation of the turbulent Prandtl number should be accounted for using 
the expression of Yakhot, et al. (1987) which gives Pr t as a function of the 
turbulent viscosity. 
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5. The effects of temperature-dependent properties should be captured by intro- 
ducing thermally expandable effects into the numerical procedure. 
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